CGEM BOSS 6.6.5.i
BESIII Offline Software System
Loading...
Searching...
No Matches
Public Member Functions | List of all members
MdcxTrackFinder Class Reference

#include <MdcxTrackFinder.h>

+ Inheritance diagram for MdcxTrackFinder:

Public Member Functions

 MdcxTrackFinder (const std::string &name, ISvcLocator *pSvcLocator)
 
virtual ~MdcxTrackFinder ()
 
StatusCode initialize ()
 
StatusCode execute ()
 
StatusCode finalize ()
 
StatusCode beginRun ()
 

Detailed Description

Definition at line 60 of file MdcxTrackFinder.h.

Constructor & Destructor Documentation

◆ MdcxTrackFinder()

MdcxTrackFinder::MdcxTrackFinder ( const std::string & name,
ISvcLocator * pSvcLocator )

Definition at line 107 of file MdcxTrackFinder.cxx.

107 :
108 Algorithm(name, pSvcLocator), m_mdcCalibFunSvc(0)
109{
110 //input
111 declareProperty("pdtFile", m_pdtFile = "pdt.table");
112 //debug control
113 declareProperty("debug", m_debug= 0);
114 declareProperty("hist", m_hist= 0);
115 declareProperty("mcHist", m_mcHist = false);
116 //cuts and control
117 declareProperty("cresol", m_cresol = 0.013);
118
119 declareProperty("getDigiFlag", m_getDigiFlag = 0);
120 declareProperty("maxMdcDigi", m_maxMdcDigi= 0);
121 declareProperty("keepBadTdc", m_keepBadTdc= 0);
122 declareProperty("dropHot", m_dropHot= 0);
123 declareProperty("keepUnmatch", m_keepUnmatch= 0);
124 declareProperty("salvageTrk", m_salvageTrk = false);
125 declareProperty("dropMultiHotInLayer",m_dropMultiHotInLayer = false);
126 declareProperty("dropTrkPt", m_dropTrkPt = -999.);
127 declareProperty("d0Cut", m_d0Cut = 999.);
128 declareProperty("z0Cut", m_z0Cut = 999.);
129
130 declareProperty("minMdcDigi", m_minMdcDigi = 0);
131 declareProperty("countPropTime",m_countPropTime = true);
132 declareProperty("addHitCut", m_addHitCut = 5.);
133 declareProperty("dropHitsSigma",m_dropHitsSigma);
134 declareProperty("helixFitCut", m_helixFitCut);
135 declareProperty("minTrkProb", m_minTrkProb = 0.01);
136 declareProperty("csmax4", m_csmax4 = 50.);
137 declareProperty("csmax3", m_csmax3 = 1.);
138 declareProperty("helixFitSigma",m_helixFitSigma= 5.);
139 declareProperty("maxRcsInAddSeg",m_maxRcsInAddSeg= 50.);
140 declareProperty("nSigAddHitTrk", m_nSigAddHitTrk= 5.);
141 declareProperty("maxProca", m_maxProca= 0.6);
142 declareProperty("doSag", m_doSag= false);
143 declareProperty("lineFit", m_lineFit = false);
144 //declareProperty("cosmicFit", m_cosmicFit= false);
145}

◆ ~MdcxTrackFinder()

MdcxTrackFinder::~MdcxTrackFinder ( )
virtual

Definition at line 150 of file MdcxTrackFinder.cxx.

150 {
151 delete m_bfield;
152}

Member Function Documentation

◆ beginRun()

StatusCode MdcxTrackFinder::beginRun ( )

Definition at line 154 of file MdcxTrackFinder.cxx.

154 {
155 // Get Mdc Detector Geometry
156 m_gm = MdcDetector::instance(m_doSag);
157 if(NULL == m_gm) return StatusCode::FAILURE;
158 MdcxHit::setMdcDetector(m_gm);
159
160 return StatusCode::SUCCESS;
161}
MdcDetector::instance
static MdcDetector * instance()
Definition MdcDetector.cxx:21
MdcxHit::setMdcDetector
static void setMdcDetector(const MdcDetector *gm)
Definition MdcxHit.cxx:68

◆ execute()

StatusCode MdcxTrackFinder::execute ( )

Definition at line 242 of file MdcxTrackFinder.cxx.

242 {
243
244
245 b_saveEvent=false;
246 setFilterPassed(b_saveEvent);
247#ifdef MDCXTIMEDEBUG
248 m_timer[0]->start();
249#endif
250 MsgStream log(msgSvc(), name());
251 log << MSG::INFO << "in execute()" << endreq;
252 StatusCode sc;
253
254 nTk = 0; t_nTdsTk = 0; //yzhang for fill
255 //------------------------------------
256 // Get event No.
257 //------------------------------------
258 SmartDataPtr<Event::EventHeader> evtHead(eventSvc(),"/Event/EventHeader");
259 if (!evtHead) {
260 log << MSG::FATAL<< "Could not retrieve event header" << endreq;
261 return StatusCode::FAILURE;
262 }
263 m_eventNo = evtHead->eventNumber();
264 //std::cout << "x evt: "<<evtHead->runNumber()<<" " << m_eventNo<< std::endl;
265 long t_evtNo = m_eventNo;
266 g_eventNo = m_eventNo;
267 //if (t_evtNo % 1000 == 0) std::cout << "x evt: " << t_evtNo << std::endl;
268 IDataManagerSvc *dataManSvc;
269 DataObject *aTrackCol;
270 DataObject *aHitCol;
271 if (!m_salvageTrk) {
272 SmartIF<IDataManagerSvc> dataManSvc(eventSvc());
273 eventSvc()->findObject("/Event/Recon/RecMdcTrackCol",aTrackCol);
274 if(aTrackCol != NULL) {
275 dataManSvc->clearSubTree("/Event/Recon/RecMdcTrackCol");
276 eventSvc()->unregisterObject("/Event/Recon/RecMdcTrackCol");
277 }
278 eventSvc()->findObject("/Event/Recon/RecMdcHitCol",aHitCol);
279 if(aHitCol != NULL) {
280 dataManSvc->clearSubTree("/Event/Recon/RecMdcHitCol");
281 eventSvc()->unregisterObject("/Event/Recon/RecMdcHitCol");
282 }
283 }
284
285 //------------------------------------
286 // Initialize track collection in TDS
287 //------------------------------------
288 DataObject *aReconEvent;
289 eventSvc()->findObject("/Event/Recon",aReconEvent);
290 if (aReconEvent==NULL) {
291 aReconEvent = new ReconEvent();
292 sc = eventSvc()->registerObject("/Event/Recon",aReconEvent);
293 if(sc != StatusCode::SUCCESS) {
294 log << MSG::FATAL << "Could not register ReconEvent" <<endreq;
295 return StatusCode::FAILURE;
296 }
297 }
298 RecMdcTrackCol* trackList;
299 eventSvc()->findObject("/Event/Recon/RecMdcTrackCol",aTrackCol);
300 if (aTrackCol) {
301 trackList = dynamic_cast<RecMdcTrackCol*> (aTrackCol);
302 }else{
303 trackList = new RecMdcTrackCol;
304 sc = eventSvc()->registerObject(EventModel::Recon::RecMdcTrackCol, trackList);
305 if(!sc.isSuccess()) {
306 log << MSG::FATAL << " Could not register RecMdcTrack collection" <<endreq;
307 return StatusCode::FAILURE;
308 }
309 }
310 RecMdcHitCol* hitList;
311 eventSvc()->findObject("/Event/Recon/RecMdcHitCol",aHitCol);
312 if (aHitCol) {
313 hitList = dynamic_cast<RecMdcHitCol*> (aHitCol);
314 }else{
315 hitList = new RecMdcHitCol;
316 sc = eventSvc()->registerObject(EventModel::Recon::RecMdcHitCol, hitList);
317 if(!sc.isSuccess()) {
318 log << MSG::FATAL << " Could not register RecMdcHit collection" <<endreq;
319 return StatusCode::FAILURE;
320 }
321 }
322
323 //------------------------------------
324 // Initialize hit collection in TDS
325 //------------------------------------
326 DataObject *pnode = 0;
327 sc = eventSvc()->retrieveObject("/Event/Hit",pnode);
328 if(!sc.isSuccess()) {
329 pnode = new DataObject;
330 sc = eventSvc()->registerObject("/Event/Hit",pnode);
331 if(!sc.isSuccess()) {
332 log << MSG::FATAL << " Could not register /Event/Hit branch " <<endreq;
333 return StatusCode::FAILURE;
334 }
335 }
336 SmartDataPtr<MdcHitCol> m_hitCol(eventSvc(), "/Event/Hit/MdcHitCol");
337 if (!m_hitCol){
338 m_hitCol= new MdcHitCol;
339 sc = eventSvc()->registerObject("/Event/Hit/MdcHitCol",m_hitCol);
340 if(!sc.isSuccess()) {
341 log << MSG::FATAL << " Could not register hit collection" <<endreq;
342 return StatusCode::FAILURE;
343 }
344 }
345
346
347 //------------------------------------
348 // Get bunch time t0 (ns) and timing
349 //------------------------------------
350 m_bunchT0 = -999.;
351 SmartDataPtr<RecEsTimeCol> aevtimeCol(eventSvc(),"/Event/Recon/RecEsTimeCol");
352 if (!aevtimeCol || aevtimeCol->size()==0) {
353 log << MSG::WARNING<< "evt "<<m_eventNo<<" Could not find RecEsTimeCol"<< endreq;
354 return StatusCode::SUCCESS;
355 }
356
357 RecEsTimeCol::iterator iter_evt = aevtimeCol->begin();
358 for(; iter_evt!=aevtimeCol->end(); iter_evt++){
359 m_bunchT0 = (*iter_evt)->getTest();
360 m_t0Stat = (*iter_evt)->getStat();
361 if ((m_t0Stat==0) || (m_bunchT0 < 0.) || (m_bunchT0 > 9999.0) ){
362 log << MSG::WARNING << "Skip evt:"<<m_eventNo<< " by t0 = "<<m_bunchT0 << endreq;
363 //return StatusCode::SUCCESS;
364 }
365 }
366 if(m_debug>1) std::cout<<name()<<" t0 "<<m_bunchT0<<" t0Stat "<<m_t0Stat<<std::endl;
367 int trigtiming=-10;
368 SmartDataPtr<TrigData> trigData(eventSvc(),"/Event/Trig/TrigData");
369 if(trigData){
370 log << MSG::INFO <<"Trigger conditions 0--43:"<<endreq;
371 for(int i = 0; i < 48; i++) {
372 log << MSG::INFO << trigData->getTrigCondName(i)<<" ---- "<<trigData->getTrigCondition(i)<< endreq;
373 }
374 for(int i = 0; i < 16; i++) log << MSG::INFO << "Trigger channel "<< i << ": " << trigData->getTrigChannel(i) << endreq;
375 m_timing=trigData->getTimingType();
376 //cout<<"-----------------trigger timing type-----------------------: "<<trigtiming<<endl;
377 log << MSG::INFO <<"Tigger Timing type: "<< trigtiming << endreq;
378 }
379
380 //------------------------------------
381 // Initialize MdcxHits
382 //------------------------------------
383 m_mdcxHits.reset();
384 uint32_t getDigiFlag = 0;
385 getDigiFlag += m_maxMdcDigi;
386 if(m_dropHot||m_salvageTrk) getDigiFlag |= MdcRawDataProvider::b_dropHot;
387 if(m_keepBadTdc) getDigiFlag |= MdcRawDataProvider::b_keepBadTdc;
388 if(m_keepUnmatch) getDigiFlag |= MdcRawDataProvider::b_keepUnmatch;
389 mdcDigiVec = m_rawDataProviderSvc->getMdcDigiVec(getDigiFlag);
390 t_nDigi = mdcDigiVec.size();
391
392 if (0 == t_nDigi ){
393 log << MSG::WARNING << " No hits in MdcDigiVec" << endreq;
394 return StatusCode::SUCCESS;
395 }
396
397
398 // fill Mc truth
399 //if(m_hist) fillMcTruth();
400
401 //skip event by hit numbe
402 if (t_nDigi < m_minMdcDigi){
403 log << MSG::WARNING << " Skip this event for MdcDigiVec.size() < "<<m_minMdcDigi << endreq;
404 return StatusCode::SUCCESS;
405 }
406 m_mdcxHits.create(mdcDigiVec, m_bunchT0, m_cresol);
407 const HepAList<MdcxHit>& dchitlist = m_mdcxHits.GetMdcxHitList();
408
409 //--------------------------------------------
410 // Make segments (MdcxSeg's) out of MdcxHit's
411 //--------------------------------------------
412 MdcxFindSegs dcsegs(dchitlist,m_debug);
413 const HepAList<MdcxSeg>& seglist = dcsegs.GetMdcxSeglist();
414 if(m_debug > 1 ){ dumpMdcxSegs(seglist);}
415 //if(m_hist){ fillMdcxSegs(seglist);}
416
417 //--------------------------------------------
418 // Make tracks (MdcxFittedHel's) out of MdcxSeg's
419 //--------------------------------------------
420 MdcxFindTracks dctrks(seglist,m_debug);
421 HepAList<MdcxFittedHel>& firsttrkl = (HepAList<MdcxFittedHel>&)dctrks.GetMdcxTrklist();
422 if(m_debug>1) dumpTrackList(firsttrkl);
423
424 //if(m_hist){ fillTrkl(firsttrkl);}
425
426 //if(m_debug>1){
427 // std::cout << "dchitlist after find tracks before MergeDups, nhits=" << dchitlist.length() << std::endl;
428 // for (int ii = 0; ii < dchitlist.length(); ii++) {
429 // dchitlist[ii]->print(std::cout, ii);
430 // }
431 // std::cout<<std::endl;
432 //}
433 MdcxMergeDups dcmergeem(firsttrkl,m_debug);
434 HepAList<MdcxFittedHel>& trkl = (HepAList<MdcxFittedHel>&)dcmergeem.GetMergedTrklist();
435
436 //if (m_debug > 1 ){
437 // cout << "MdcxTrackFinder: after MergeDups, have "
438 // << trkl.length() << " track(s). nhits=" << dchitlist.length() << endl;
439 // for (int ii = 0; ii < dchitlist.length(); ii++) {
440 // dchitlist[ii]->print(std::cout, ii);
441 // }
442 // std::cout<<std::endl;
443 //}
444
445 //---------------------------------------------------------
446 // Put my tracks into official fitter and store to TDS
447 //----------------------------------------------------------
448
449
450 sc = FitMdcxTrack(trkl, dchitlist, m_hitCol, trackList, hitList);
451 if (!sc.isSuccess()) {return StatusCode::SUCCESS;}
452 t_nTdsTk = trackList->size();
453
454 t_nTkTot += trackList->size();
455 if(t_nTdsTk<20) t_nTkNum[t_nTdsTk]++;
456
457#ifdef MDCXTIMEDEBUG
458 m_timer[0]->stop();
459 m_xt4time = m_timer[0]->elapsed();
460 m_xtupleEvt->write();
461#endif
462
463 if(m_hist) fillEvent();
464 if (m_debug > 0) {
465 DataObject* pNode;
466 eventSvc()->retrieveObject("/Event/Recon/RecMdcTrackCol",pNode);
467 RecMdcTrackCol *tmpTrackCol = dynamic_cast<RecMdcTrackCol*> (pNode);
468 eventSvc()->retrieveObject("/Event/Recon/RecMdcHitCol",pNode);
469 int nTdsTk = 0;
470 if(tmpTrackCol) nTdsTk = tmpTrackCol->size();
471
472 //if (t_evtNo % 1000 == 0) {
473 std::cout<< "MdcxTrackFinder: evtNo "<< m_eventNo << " t0="<<m_bunchT0
474 <<" Found " <<trkl.length()
475 <<" keep "<< t_nTdsTk
476 <<" finialy keep "<< nTdsTk;
477
478 int ndelete =0; trkl.length() - trackList->size();
479 if( ndelete>0 ) std::cout <<" delete "<< ndelete;
480 std::cout <<" track(s)" <<endl;
481 //}
482
483 if(m_debug>1)dumpTdsTrack(tmpTrackCol);
484 if(m_debug>1)dumpTrack(tmpTrackCol);
485 //dumpTdsHits(tmpHitCol);
486
487 }
488 if((trackList->size()!=4) ) b_saveEvent = true;
489 setFilterPassed(b_saveEvent);
490 return StatusCode::SUCCESS;
491}
MdcHitCol
ObjectVector< MdcHit > MdcHitCol
Definition MdcHit.h:129
m_xt4time
NTuple::Item< double > m_xt4time
Definition MdcxHistItem.h:122
m_xtupleEvt
NTuple::Tuple * m_xtupleEvt
Definition MdcxHistItem.h:64
RecMdcHitCol
ObjectVector< RecMdcHit > RecMdcHitCol
Definition RecMdcHit.h:99
RecMdcTrackCol
ObjectVector< RecMdcTrack > RecMdcTrackCol
Definition RecMdcTrack.h:102
msgSvc
IMessageSvc * msgSvc()
Definition ServiceAccessor.h:13
MdcxHits::GetMdcxHitList
const HepAList< MdcxHit > & GetMdcxHitList()
Definition MdcxHits.h:40
MdcxHits::reset
void reset()
Definition MdcxHits.cxx:55
MdcxHits::create
void create(MdcDigiVec digiVec, float c0=0.0, float cresol=0.0180)
Definition MdcxHits.cxx:59
RawDataProviderSvc::getMdcDigiVec
MdcDigiVec & getMdcDigiVec(uint32_t control=0)
Definition RawDataProviderSvc.h:29
EventModel::Recon::RecMdcTrackCol
_EXTERN_ std::string RecMdcTrackCol
Definition EventModel.h:92
EventModel::Recon::RecMdcHitCol
_EXTERN_ std::string RecMdcHitCol
Definition EventModel.h:91
g_eventNo
NTuple::Item< long > g_eventNo
Definition ntupleItem.h:22

◆ finalize()

StatusCode MdcxTrackFinder::finalize ( )

Definition at line 493 of file MdcxTrackFinder.cxx.

493 {
494 MsgStream log(msgSvc(), name());
495 log << MSG::INFO << "in finalize()" << endreq;
496
497 std::cout<< " --after " << name() << " keep " << t_nTkTot<< " tracks "<<std::endl;
498 for(int i=0; i<20; i++){
499 if(t_nTkNum[i]>0)std::cout<< " nTk="<< i << " "<< t_nTkNum[i]<< std::endl;
500 }
501
502 //tot evtNo, trkNum
503 return StatusCode::SUCCESS;
504}

◆ initialize()

StatusCode MdcxTrackFinder::initialize ( )

Definition at line 165 of file MdcxTrackFinder.cxx.

165 {
166 MsgStream log(msgSvc(), name());
167 log << MSG::INFO << "in initialize()" << endreq;
168
169 t_nTkTot = 0;
170 for (int i=0;i<20;i++) t_nTkNum[i]=0;
171
172 //m_flags.readPar(m_paramFile);
173#ifdef MDCXTIMEDEBUG
174 StatusCode tsc = service( "BesTimerSvc", m_timersvc);
175 if( tsc.isFailure() ) {
176 log << MSG::WARNING << name() << ": Unable to locate BesTimer Service" << endreq;
177 return StatusCode::FAILURE;
178 }
179 m_timer[0] = m_timersvc->addItem("Execution");
180 m_timer[0]->propName("Execution");
181#endif
182
183 if (m_helixFitCut.size() == 43){
184 for(int i=0; i<43; i++){
185 //MdcTrkReconCut_helix_fit[i] = m_helixFitCut[i];
186 TrkHelixFitter::nSigmaCut[i] = m_helixFitCut[i];
187 }
188 }
189 MdcxParameters::debug = m_debug;
190 MdcxParameters::minTrkProb = m_minTrkProb;
191 MdcxParameters::csmax4 = m_csmax4;
192 MdcxParameters::csmax3 = m_csmax3;
193 MdcxParameters::helixFitSigma = m_helixFitSigma;
194 MdcxParameters::maxRcsInAddSeg= m_maxRcsInAddSeg;
195 MdcxParameters::nSigAddHitTrk = m_nSigAddHitTrk;
196 MdcxParameters::maxProca = m_maxProca;
197 TrkHelixFitter::m_debug = (m_debug>7);
198 Pdt::readMCppTable(m_pdtFile);
199 MdcxFittedHel::debug = m_debug;
200
201
202 // Get MdcCalibFunSvc
203 IMdcCalibFunSvc* imdcCalibSvc;
204 StatusCode sc = service ("MdcCalibFunSvc", imdcCalibSvc);
205 m_mdcCalibFunSvc = dynamic_cast<MdcCalibFunSvc*>(imdcCalibSvc);
206 if ( sc.isFailure() ){
207 log << MSG::FATAL << "Could not load MdcCalibFunSvc!" << endreq;
208 return StatusCode::FAILURE;
209 }
210 MdcxHit::setMdcCalibFunSvc(m_mdcCalibFunSvc);
211 MdcxHit::setCountPropTime(m_countPropTime);
212
213 // Get RawDataProviderSvc
214 IRawDataProviderSvc* iRawDataProvider;
215 sc = service ("RawDataProviderSvc", iRawDataProvider);
216 if ( sc.isFailure() ){
217 log << MSG::FATAL << "Could not load RawDataProviderSvc!" << endreq;
218 return StatusCode::FAILURE;
219 }
220 m_rawDataProviderSvc = dynamic_cast<RawDataProviderSvc*>(iRawDataProvider);
221
222
223 //Initailize magnetic filed
224 sc = service ("MagneticFieldSvc",m_pIMF);
225 if(sc!=StatusCode::SUCCESS) {
226 log << MSG::ERROR << "Unable to open Magnetic field service"<<endreq;
227 }
228 m_bfield = new BField(m_pIMF);
229 log << MSG::INFO << "field z = "<<m_bfield->bFieldNominal()<< endreq;
230 m_context = new TrkContextEv(m_bfield);
231
232 if (m_hist) {bookNTuple();}
233 if (m_dropHitsSigma.size()==43){
234 for (int ii=0;ii<43;ii++) {
235 MdcxParameters::dropHitsSigma[ii]=m_dropHitsSigma[ii];
236 }
237 }
238
239 return StatusCode::SUCCESS;
240}
MdcxFittedHel::debug
static int debug
Definition MdcxFittedHel.h:116
MdcxHit::setMdcCalibFunSvc
static void setMdcCalibFunSvc(const MdcCalibFunSvc *calibSvc)
Definition MdcxHit.cxx:60
MdcxHit::setCountPropTime
static void setCountPropTime(bool countPropTime)
Definition MdcxHit.cxx:64
MdcxParameters::maxProca
static double maxProca
Definition MdcxParameters.h:14
MdcxParameters::debug
static int debug
Definition MdcxParameters.h:44
MdcxParameters::maxRcsInAddSeg
static double maxRcsInAddSeg
Definition MdcxParameters.h:21
MdcxParameters::minTrkProb
static double minTrkProb
Definition MdcxParameters.h:15
MdcxParameters::dropHitsSigma
static float dropHitsSigma[43]
Definition MdcxParameters.h:35
MdcxParameters::csmax4
static double csmax4
Definition MdcxParameters.h:29
MdcxParameters::helixFitSigma
static double helixFitSigma
Definition MdcxParameters.h:26
MdcxParameters::csmax3
static double csmax3
Definition MdcxParameters.h:30
MdcxParameters::nSigAddHitTrk
static double nSigAddHitTrk
Definition MdcxParameters.h:22
Pdt::readMCppTable
static void readMCppTable(std::string filenm)
Definition Pdt.cxx:291
TrkHelixFitter::m_debug
static bool m_debug
Definition TrkHelixFitter.h:40
TrkHelixFitter::nSigmaCut
static double nSigmaCut[43]
Definition TrkHelixFitter.h:37
The documentation for this class was generated from the following files: