BOSS 7.0.5
BESIII Offline Software System
Loading...
Searching...
No Matches
MdcTrkRecon.cxx
Go to the documentation of this file.
1#include "MdcTrkRecon/MdcTrkRecon.h"
2#include "GaudiKernel/MsgStream.h"
3#include "GaudiKernel/AlgFactory.h"
4#include "GaudiKernel/ISvcLocator.h"
5#include "GaudiKernel/DataSvc.h"
6#include "GaudiKernel/SmartDataPtr.h"
7#include "GaudiKernel/IDataProviderSvc.h"
8#include "GaudiKernel/IDataManagerSvc.h"
9#include "GaudiKernel/PropertyMgr.h"
10#include "GaudiKernel/Bootstrap.h"
11#include "GaudiKernel/INTupleSvc.h"
12#include "EventModel/EventHeader.h"
13#include "ReconEvent/ReconEvent.h"
14#include "MdcRecEvent/RecMdcHit.h"
15#include "MdcRecEvent/RecMdcTrack.h"
16#include "RawEvent/RawDataUtil.h"
17#include "RawDataProviderSvc/IRawDataProviderSvc.h"
18#include "RawDataProviderSvc/MdcRawDataProvider.h"
19#include "Identifier/MdcID.h"
20#include "Identifier/Identifier.h"
21#include "MdcRawEvent/MdcDigi.h"
22#include "MdcData/MdcHit.h"
23#include "MdcData/MdcHitMap.h"
24#include "MdcData/MdcHitUse.h"
25#include "MdcData/MdcHitOnTrack.h"
26#include "MdcData/MdcRecoHitOnTrack.h"
27#include "MdcGeom/MdcDetector.h"
28#include "EvTimeEvent/RecEsTime.h"
29
30#include "MdcTrkRecon/MdcSegData.h"
31#include "MdcTrkRecon/MdcSegList.h"
32#include "MdcTrkRecon/MdcSegFinder.h"
33#include "MdcTrkRecon/MdcTrackList.h"
34#include "MdcTrkRecon/MdcTrackListCsmc.h"
35#include "TrkFitter/TrkContextEv.h"
36#include "BField/BField.h"
37#include "MdcTrkRecon/MdcHistItem.h"
38#include "MdcTrkRecon/MdcSeg.h"
39#include "MdcTrkRecon/GmsList.h"
40#include "MdcGeom/Constants.h"
41#include "MdcRecoUtil/Pdt.h"
42#include "MdcTrkRecon/MdcTrack.h"
43#include "TrkBase/TrkExchangePar.h"
44#include "TrkBase/TrkRep.h"
45#include "TrkBase/TrkDifTraj.h"
46#include "TrkFitter/TrkHelixFitter.h"
47
48#include "McTruth/McParticle.h"
49#include "McTruth/MdcMcHit.h"
50#include "MdcPrintSvc/IMdcPrintSvc.h"
51#include "MdcTrkRecon/MdcSegUsage.h"
52
53
54#ifdef MDCPATREC_TIMETEST
55#include "BesTimerSvc/IBesTimerSvc.h"
56#include "BesTimerSvc/BesTimerSvc.h"
57#include "BesTimerSvc/BesTimer.h"
58#endif
59
60#include <vector>
61#include <iostream>
62#include <fstream>
63
64using namespace std;
65using namespace Event;
66
67//-----------------------------------------
68
69/////////////////////////////////////////////////////////////////////////////
70
71MdcTrkRecon::MdcTrkRecon(const std::string& name, ISvcLocator* pSvcLocator) :
72 Algorithm(name, pSvcLocator),
73 m_hitData(0), m_segs(0), m_tracks(0), m_segFinder(0)
74{
75 //Declare the properties--------------------------------
76 declareProperty("FittingMethod", m_fittingMethod = 2);
77 declareProperty("ConfigFile", m_configFile = "MDCConfig.xml");
78 declareProperty("OnlyUnusedHits", m_onlyUnusedHits = 0);
79 declareProperty("PoisonHits", m_poisonHits = 0);
80 declareProperty("doLineFit", m_doLineFit = false);
81 declareProperty("paramFile", m_paramFile = "params");
82 declareProperty("pdtFile", m_pdtFile = "pdt.table");
83 declareProperty("allHit", m_allHit = 1);
84 declareProperty("hist", m_hist= false);
85 declareProperty("recForEsTime", m_recForEsTime= false);
86 declareProperty("tryBunch", m_tryBunch = false);
87 declareProperty("d0Cut", m_d0Cut= -999.);
88 declareProperty("z0Cut", m_z0Cut= -999.);
89 declareProperty("dropTrkPt", m_dropTrkPt = -999.);
90 declareProperty("debug", m_debug= 0);
91 declareProperty("mcHist", m_mcHist= false);
92 declareProperty("fieldCosmic", m_fieldCosmic= false);
93 declareProperty("doSag", m_doSagFlag= false);
94 declareProperty("arbitrateHits", m_arbitrateHits = true);
95
96 declareProperty("helixHitsSigma", m_helixHitsSigma);
97 //for fill hist
98 declareProperty("getDigiFlag", m_getDigiFlag = 0);
99 declareProperty("maxMdcDigi", m_maxMdcDigi = 0);
100 declareProperty("keepBadTdc", m_keepBadTdc = 0);
101 declareProperty("dropHot", m_dropHot= 0);
102 declareProperty("keepUnmatch", m_keepUnmatch = 0);
103 declareProperty("minMdcDigi", m_minMdcDigi = 0);
104 declareProperty("selEvtNo", m_selEvtNo);
105 declareProperty("combineTracking",m_combineTracking = false);//yzhang 2010-05-28
106 declareProperty("noInner",m_noInner = false);//yzhang 2016-10-12
107
108#ifdef MDCPATREC_RESLAYER
109 declareProperty("resLayer",m_resLayer= -1);
110#endif
111
112}
113
114MdcTrkRecon::~MdcTrkRecon()
115{
116 m_segs.reset(0);
117 m_segFinder.reset(0);
118 m_hitData.reset(0);
119 m_tracks.reset(0);
120}
121
122// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
123StatusCode MdcTrkRecon::beginRun(){
124 MsgStream log(msgSvc(), name());
125 log << MSG::INFO << "in beginRun()" << endreq;
126
127 //Detector geometry
128 _gm = MdcDetector::instance(m_doSagFlag);
129
130 if(NULL == _gm) return StatusCode::FAILURE;
131
132 //Initialize segList
133 m_segs.reset( new MdcSegList(_gm->nSuper(), m_flags.segPar) );
134
135 return StatusCode::SUCCESS;
136}
137
138// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
139StatusCode MdcTrkRecon::initialize(){
140
141 MsgStream log(msgSvc(), name());
142 log << MSG::INFO << "in initialize()" << endreq;
143 StatusCode sc;
144
145 //Pdt
146 Pdt::readMCppTable(m_pdtFile);
147
148 //Flag and Pars
149 m_flags.readPar(m_paramFile);
150 m_flags.lHist = m_hist;
151 m_flags.setDebug(m_debug);
152 for(int i=0;i<43;i++) TrkHelixFitter::nSigmaCut[i] = m_helixHitsSigma[i];
153 if (!m_doLineFit) MdcTrackListBase::m_d0Cut = m_d0Cut;
154 if (!m_doLineFit) MdcTrackListBase::m_z0Cut = m_z0Cut;
155 MdcTrackListBase::m_ptCut = m_dropTrkPt;
156 if (m_debug>0) {
157 std::cout<< "MdcTrkRecon d0Cut "<<m_d0Cut<< " z0Cut "<<
158 m_z0Cut<<" ptCut "<<m_dropTrkPt << std::endl;
159 }
160
161 //Initailize magnetic filed
162 sc = service ("MagneticFieldSvc",m_pIMF);
163 if(sc!=StatusCode::SUCCESS) {
164 log << MSG::ERROR << name() << " Unable to open magnetic field service"<<endreq;
165 }
166 m_bfield = new BField(m_pIMF);
167 log << MSG::INFO <<name() << " field z = "<<m_bfield->bFieldNominal()<< endreq;
168
169 //Initialize RawDataProviderSvc
170 IRawDataProviderSvc* irawDataProviderSvc;
171 sc = service ("RawDataProviderSvc", irawDataProviderSvc);
172 m_rawDataProviderSvc = dynamic_cast<RawDataProviderSvc*> (irawDataProviderSvc);
173 if ( sc.isFailure() ){
174 log << MSG::FATAL << name()<<" Could not load RawDataProviderSvc!" << endreq;
175 return StatusCode::FAILURE;
176 }
177
178 //Initailize MdcPrintSvc
179 IMdcPrintSvc* imdcPrintSvc;
180 sc = service ("MdcPrintSvc", imdcPrintSvc);
181 m_mdcPrintSvc = dynamic_cast<MdcPrintSvc*> (imdcPrintSvc);
182 if ( sc.isFailure() ){
183 log << MSG::FATAL << "Could not load MdcPrintSvc!" << endreq;
184 return StatusCode::FAILURE;
185 }
186
187 //Initialize hitdata
188 m_hitData.reset( new MdcSegData( ignoringUsedHits() ));
189
190 //Initialize segFinder
191 if (m_flags.findSegs) {
192 m_segFinder.reset( new MdcSegFinder(m_flags.segPar.useAllAmbig) );
193 }
194
195 //Initialize trackList
196 if ( m_doLineFit ){
197 m_tracks.reset( new MdcTrackListCsmc(m_flags.tkParTight) );
198 } else {
199 m_tracks.reset( new MdcTrackList(m_flags.tkParTight) );
200 }
201 m_tracks->setNoInner(m_noInner);
202
203 //Book NTuple and Histograms
204 if (m_flags.lHist){
205 StatusCode sc = bookNTuple();
206 if (!sc.isSuccess()) { m_flags.lHist = 0;}
207 }
208
209 //yzhang for HelixFitter debug
210 if(7 == m_flags.debugFlag())TrkHelixFitter::m_debug = true;
211
212 t_iExexute = 0;
213 return StatusCode::SUCCESS;
214}
215
216// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
217StatusCode MdcTrkRecon::execute() {
218 setFilterPassed(false);
219
220 MsgStream log(msgSvc(), name());
221 log << MSG::INFO << "in execute()" << endreq;
222
223 //Initialize
224 if(m_hist>0){
225 for (int ii=0;ii<43;ii++){
226 for (int jj=0;jj<288;jj++){
227 haveDigiTk[ii][jj] = -2;
228 haveDigiPt[ii][jj] = -2;
229 haveDigiTheta[ii][jj] = -999.;
230 haveDigiPhi[ii][jj] = -999.;
231 haveDigiDrift[ii][jj] = -999.;
232 haveDigiAmbig[ii][jj] = -999.;
233 recHitMap[ii][jj] = 0;
234 mcDrift[ii][jj] = -99.;
235 mcX[ii][jj] = -99.;
236 mcY[ii][jj] = -99.;
237 mcZ[ii][jj] = -99.;
238 mcLR[ii][jj] = -99.;
239 hitPoisoned[ii][jj]=-99;
240 }
241 }
242 }
243
244
245 TrkContextEv context(m_bfield);
246#ifdef MDCPATREC_TIMETEST
247 m_timer[1]->start();
248 ti_nHit=0;
249#endif
250 //------------------read event header--------------
251 SmartDataPtr<Event::EventHeader> evtHead(eventSvc(),"/Event/EventHeader");
252 if (!evtHead) {
253 log << MSG::FATAL<< "Could not retrieve event header" << endreq;
254 return( StatusCode::FAILURE);
255 }
256 t_eventNo = evtHead->eventNumber();
257 if(m_debug!=0)std::cout<<t_iExexute<<" run "<<evtHead->runNumber()<<" evt "<<t_eventNo<<std::endl;
258 t_iExexute++;
259
260 if (m_selEvtNo.size() >0){
261 std::vector<int>::iterator it = m_selEvtNo.begin();
262 for (; it!=m_selEvtNo.end(); it++){
263 if((*it)== t_eventNo) setFilterPassed(true);
264 }
265 }
266 //------------------get event start time-----------
267 double t0 = 0.;
268 t_t0 = -1.;
269 t_t0Stat = -1;
270 bool iterateT0 = false;
271 const int nBunch = 3;//3 bunch/event
272 double iBunch = 0 ; //form 0 to 2
273 double BeamDeltaTime = 24. / nBunch;// 8ns
274 SmartDataPtr<RecEsTimeCol> evTimeCol(eventSvc(),"/Event/Recon/RecEsTimeCol");
275 if (!evTimeCol || evTimeCol->size()==0) {
276 log << MSG::WARNING <<" run "<<evtHead->runNumber() <<" evt "<<t_eventNo<< " Could not find RecEsTimeCol" << endreq;
277 if(m_fieldCosmic){//yzhang temp for bes3 csmc test
278 return StatusCode::SUCCESS;
279 }
280 }
281 RecEsTimeCol::iterator iter_evt = evTimeCol->begin();
282 // skip RecEsTimeCol no rec data
283 if (iter_evt != evTimeCol->end()){
284 t_t0Stat = (*iter_evt)->getStat();
285 t_t0 = (*iter_evt)->getTest();
286
287 if (m_doLineFit){
288 //t0 -= m_t0Const;
289 if (abs(t_t0)<0.00001 || (t_t0 < 0) || (t_t0 > 9999.0) ){
290 log << MSG::WARNING << "Skip with t0 = "<<t_t0 << endreq;
291 return StatusCode::SUCCESS;//for bes3 cosmic test
292 }
293 }
294 t0 = t_t0*1.e-9;
295 if(m_debug>0) std::cout<<name()<<" t0 "<<t0<<" "<<std::endl;
296 }
297
298
299restart:
300 if ( !m_recForEsTime && m_tryBunch) {
301 if ( t_t0Stat < 10 ) return StatusCode::SUCCESS;
302 }
303 if ( m_tryBunch ){
304 if ( iBunch > (nBunch - 1) ) return StatusCode::SUCCESS;
305 if ( t_t0Stat < 0 ){ iterateT0 = true; }
306 if ( iterateT0 ){
307 //double EventDeltaTime = 24.;//24ns/event
308 if ( (t_t0Stat > -1) && (fabs( iBunch * BeamDeltaTime - t_t0) < 2.) ){
309 ++iBunch;
310 goto restart;
311 }else{ t0 = iBunch * BeamDeltaTime *1.e-9;}
312 ++iBunch;
313 }
314 }
315
316 SmartDataPtr<MdcHitMap> hitMap(eventSvc(),"/Event/Hit/MdcHitMap");
317 if (!hitMap) {
318 log << MSG::FATAL << "Could not retrieve MDC hit map" << endreq;
319 return( StatusCode::FAILURE);
320 }
321 //---------------------Read an event--------------
322 SmartDataPtr<MdcHitCol> hitCol(eventSvc(),"/Event/Hit/MdcHitCol");
323 if (!hitCol) {
324 log << MSG::FATAL << "Could not retrieve MDC hit list" << endreq;
325 return( StatusCode::FAILURE);
326 }
327 StatusCode sc;
328
329 //---------- register RecMdcTrackCol and RecMdcHitCol to tds---------
330 DataObject *aReconEvent;
331 eventSvc()->findObject("/Event/Recon",aReconEvent);
332 if(aReconEvent==NULL) {
333 aReconEvent = new ReconEvent();
334 StatusCode sc = eventSvc()->registerObject("/Event/Recon",aReconEvent);
335 if(sc!=StatusCode::SUCCESS) {
336 log << MSG::FATAL << "Could not register ReconEvent" <<endreq;
337 return( StatusCode::FAILURE);
338 }
339 }
340
341 DataObject *aTrackCol;
342 DataObject *aRecHitCol;
343 //yzhang 2010-05-28
344 SmartIF<IDataManagerSvc> dataManSvc(eventSvc());
345 if(!m_combineTracking){
346 eventSvc()->findObject("/Event/Recon/RecMdcTrackCol",aTrackCol);
347 if(aTrackCol != NULL) {
348 dataManSvc->clearSubTree("/Event/Recon/RecMdcTrackCol");
349 eventSvc()->unregisterObject("/Event/Recon/RecMdcTrackCol");
350 }
351 eventSvc()->findObject("/Event/Recon/RecMdcHitCol",aRecHitCol);
352 if(aRecHitCol != NULL) {
353 dataManSvc->clearSubTree("/Event/Recon/RecMdcHitCol");
354 eventSvc()->unregisterObject("/Event/Recon/RecMdcHitCol");
355 }
356 }
357
358 RecMdcTrackCol* trackList;
359 eventSvc()->findObject("/Event/Recon/RecMdcTrackCol",aTrackCol);
360 if (aTrackCol) {
361 trackList = dynamic_cast<RecMdcTrackCol*> (aTrackCol);
362 }else{
363 trackList = new RecMdcTrackCol;
364 sc = eventSvc()->registerObject(EventModel::Recon::RecMdcTrackCol, trackList);
365 if(!sc.isSuccess()) {
366 log << MSG::FATAL << " Could not register RecMdcTrack collection" <<endreq;
367 return StatusCode::FAILURE;
368 }
369 }
370 RecMdcHitCol* hitList;
371 eventSvc()->findObject("/Event/Recon/RecMdcHitCol",aRecHitCol);
372 if (aRecHitCol) {
373 hitList = dynamic_cast<RecMdcHitCol*> (aRecHitCol);
374 }else{
375 hitList = new RecMdcHitCol;
376 sc = eventSvc()->registerObject(EventModel::Recon::RecMdcHitCol, hitList);
377 if(!sc.isSuccess()) {
378 log << MSG::FATAL << " Could not register RecMdcHit collection" <<endreq;
379 return StatusCode::FAILURE;
380 }
381 }
382
383 if(m_debug>0) std::cout<<name()<<" t0 "<<t0<<" "<<std::endl;
384 m_hitData->loadevent(hitCol, hitMap, t0);
385 t_nDigi = hitCol->size();
386
387 if (poisoningHits()) {
388 m_hitData->poisonHits(_gm,m_debug);
389 }
390
391 if((m_hist>0) && m_tupleWireEff){
392 for(int i=0;i<m_hitData->nhits();i++){
393 const MdcHit* thisHit = m_hitData->hit(i);
394 int l = thisHit->layernumber();
395 int w = thisHit->wirenumber();
396 if(m_hitData->segUsage().get(thisHit)->deadHit()){
397 hitPoisoned[l][w] = 1;
398 }else{
399 hitPoisoned[l][w] = 0;
400 }
401 }
402 }
403
404#ifdef MDCPATREC_TIMETEST
405 SmartDataPtr<Event::McParticleCol> mcParticleCol(eventSvc(),"/Event/MC/McParticleCol");
406 if(!mcParticleCol){
407 log << MSG::INFO << "Could not retrieve McParticelCol" << endreq;
408 }
409 m_mcTkNum = mcParticleCol->size();
410#endif
411
412
413 if(m_debug>1) dumpDigi();
414 //--------------------------------------------------------------------------
415 // Segment finding
416 //--------------------------------------------------------------------------
417 int nSegs = 0;
418 if (m_flags.findSegs) {
419 // Form track segments in superlayers
420 nSegs = m_segFinder->createSegs(_gm, *m_segs, m_hitData->segUsage(),
421 m_hitData->hitMap(), t0);
422 }
423 log << MSG::INFO << " Found " << nSegs << " segments" << endreq;
424 if (m_debug>1){
425 std::cout<<"------------------------------------------------"<<std::endl;
426 std::cout<<"==event "<<t_eventNo<< " Found " << nSegs << " segments" << std::endl;
427 m_segs->plot();
428 }
429
430 if (m_flags.lHist||m_flags.segPar.lPrint) {
431 fillSegList();
432 }
433
434 //--------------------------------------------------------------------------
435 // Combine segments to form track
436 //--------------------------------------------------------------------------
437 int nTracks = 0;
438 int nDeleted = 0;
439 if (m_flags.findTracks && m_flags.findSegs) {
440 if (nSegs > 2) {
441 nTracks = m_tracks->createFromSegs(m_segs.get(), m_hitData->hitMap(),
442 _gm, context, t0);
443 }
444
445 if(m_arbitrateHits) nDeleted = m_tracks->arbitrateHits();
446 int nKeep = nTracks - nDeleted;
447
448 if (m_debug>0 && (nTracks - nDeleted)>0){
449 std::cout<< "MdcTrkRecon: evt "<< setw(5)<<t_eventNo
450 <<" nDigi "<<setw(4)<<t_nDigi<< " t0 "<<setw(5)<<t_t0
451 << " keep "<< nKeep <<" track(s)";
452 if (nDeleted!=0){ std::cout <<",deleted " <<nDeleted; }
453 std::cout<< std::endl;//yzhang debug
454 }else{
455 if (m_debug>0){
456 std::cout<< "MdcTrkRecon: evt "<< setw(5)<<t_eventNo
457 <<" nDigi "<<setw(4)<<t_nDigi<< " t0 "<<setw(5)<<t_t0
458 <<" found 0 track "<< endl;
459 }
460 }
461
462 if (m_flags.lHist) t_nRecTk = nKeep;
463
464#ifdef MDCPATREC_RESLAYER
465 m_tracks->setResLayer(m_resLayer);
466#endif
467 if (m_flags.lHist){ fillTrackList(); }
468 // Stick the found tracks into the list in TDS
469 m_tracks->store(trackList, hitList);
470 //if (m_debug>1) { dumpTdsTrack(trackList); }
471 if(m_debug>1){
472 std::cout<<name()<<" print track list "<<std::endl;
473 m_mdcPrintSvc->printRecMdcTrackCol();
474 }
475
476 //if(nKeep != (int)trackList->size()) std::cout<<__FILE__<<" "<<__LINE__<<"!!!!!!!!!!!!!!!!! nKeep != nTdsTk"<< std::endl;
477#ifdef MDCPATREC_TIMETEST
478 RecMdcTrackCol::iterator iter = trackList->begin();
479 for (;iter != trackList->end(); iter++) {
480 MdcTrack* tk = *iter;
481 ti_nHit += tk->getNhits();
482 }
483#endif
484
485 }
486 //-------------End track-finding------------------
487 m_segs->destroySegs();
488
489 //if ( t_eventNo% 1000 == 0) {
490 //std::cout<<"evt "<<t_eventNo<< " Found " << trackList->size() << " track(s)" <<endl;//yzhang debug
491 //}
492 if(m_flags.lHist && m_mcHist) fillMcTruth();
493#ifdef MDCPATREC_TIMETEST
494 m_timer[1]->stop();
495 //cout << "m_timer[1]->elapsed()::" << m_timer[1]->elapsed() << endl;
496 //cout << "m_timer[1]->mean()::" << m_timer[1]->mean() << endl;
497 m_eventTime = m_timer[1]->elapsed();
498 m_t5recTkNum = m_tracks->length();
499 m_t5EvtNo = t_eventNo;
500 m_t5nHit = ti_nHit;
501 m_t5nDigi = t_nDigi;
502 sc = m_tupleTime->write();
503#endif
504 // for event start time, if track not found try another t0
505 if ( m_tryBunch ){
506 if ( nTracks <1 ){ iterateT0 = true;
507 goto restart;
508 }
509 }
510 if (m_flags.lHist ) fillEvent();
511
512 return StatusCode::SUCCESS;
513}
514
515// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
516StatusCode MdcTrkRecon::finalize() {
517 MsgStream log(msgSvc(), name());
518 log << MSG::INFO << "in finalize()" << endreq;
519
520 m_tracks.reset(0);
521#ifdef MDCPATREC_TIMETEST
522 m_timersvc->print();
523#endif
524 return StatusCode::SUCCESS;
525}
526
527StatusCode MdcTrkRecon::bookNTuple(){
528 MsgStream log(msgSvc(), name());
529 StatusCode sc = StatusCode::SUCCESS;
530
531
532 //--------------book histogram------------------
533 h_sfHit = histoSvc()->book( "sfHit", "Hits after segment finding",46,-1.,44. );
534 //g_actHit = histoSvc()->book( "actHit", "Active hits",46,-1.,44. );
535 //g_hitEff = histoSvc()->book( "hitEff", "Hit efficiency in the event",100,-0.5,1.2 );
536 g_residVsLayer = histoSvc()->book( "residVsLayer", "Residual vs. layer",60,-5,50.,100,-0.5,1.2 );
537 g_residVsDoca = histoSvc()->book( "residVsDoca", "Residial vs. doca",50,-2.,2.,100,-0.5,1.2 );
538 //segment
539 //g_segChi2 = histoSvc()->book( "segChi2", "chi2 per dof in segment finding",101,-1.,100. );
540 g_cirTkChi2 = histoSvc()->book( "cirTkChi2", "chi2 per dof in circle finding",21,-1. ,20. );
541 g_3dTkChi2 = histoSvc()->book( "chi2Helix", "maxChisq dof helix fit" ,51.,-1.,50);
542 g_nSigAdd = histoSvc()->book( "nSigAdd", "max allowed n sigma to add hit to existing segment",50,0,50 );
543 g_z0Cut = histoSvc()->book( "z0Cut", "z0 for including stereo segs in cluster",100,0,600 );
544 g_ctCut = histoSvc()->book( "ctCut", "ct for including stereo segs in cluster",100,-50,50 );
545 g_delCt = histoSvc()->book( "delCt", "del ct cut forming seg group",100,-20,20 );
546 g_delZ0 = histoSvc()->book( "delZ0", "del z0 cut forming seg group",100,-60,60 );
547 g_phiDiff = histoSvc()->book( "phiDiff", "phidiff mult for drop dup seg",100,-0.5,6.5 );
548 g_slopeDiff = histoSvc()->book( "slopeDiff", "slopediff for drop dup seg",100,-0.3,0.3 );
549 g_maxSegChisqAx = histoSvc()->book( "maxSegChisqAx", "max chisqDof ax seg combine" ,1000,0.,1000.);
550 g_maxSegChisqSt = histoSvc()->book( "maxSegChisqSt", "max chisqDof st seg combine" ,200,-10.,200.);
551 g_pickHitWire= histoSvc()->book( "pickHitWire", "nWire of pickHit" ,600,-288,288 );
552
553 //for(int i=0;i<11;i++){
554 //char title1[80],title2[80];
555 //sprintf(title1,"segChi2Pat3%d",i);
556 //sprintf(title2,"segChi2Pat4%d",i);
557 //g_segChi2SlayPat3[i] = histoSvc()->book( title1, "chi2/dof per layer of pat3",710,-10.,700. );
558 //g_segChi2SlayPat4[i] = histoSvc()->book( title2, "chi2/dof per layer of pat4",710,-10.,700. );
559 //}
560
561 //book tuple of wire efficiency
562 NTuplePtr ntWireEff(ntupleSvc(), "MdcWire/mc");
563 if ( ntWireEff ) { m_tupleWireEff = ntWireEff;}
564 else {
565 m_tupleWireEff = ntupleSvc()->book ("MdcWire/mc", CLID_ColumnWiseTuple, "MdcWire");
566 if ( m_tupleWireEff ) {
567 sc = m_tupleWireEff->addItem ("tkId", m_we_tkId);
568 sc = m_tupleWireEff->addItem ("pdg", m_we_pdg);
569 sc = m_tupleWireEff->addItem ("primary", m_we_primary);
570 sc = m_tupleWireEff->addItem ("layer", m_we_layer);
571 sc = m_tupleWireEff->addItem ("wire", m_we_wire);
572 sc = m_tupleWireEff->addItem ("gwire", m_we_gwire);
573 sc = m_tupleWireEff->addItem ("poisoned", m_we_poisoned);
574 sc = m_tupleWireEff->addItem ("seg", m_we_seg);
575 sc = m_tupleWireEff->addItem ("track", m_we_track);
576 sc = m_tupleWireEff->addItem ("pt", m_we_pt);
577 sc = m_tupleWireEff->addItem ("theta", m_we_theta);
578 sc = m_tupleWireEff->addItem ("phi", m_we_phi);
579 //sc = m_tupleWireEff->addItem ("tdc", m_we_tdc);
580 } else {
581 log << MSG::ERROR << "Cannot book tuple MdcWire/mc" << endmsg;
582 }
583 }
584
585 //book tuple of test
586 NTuplePtr ntt(ntupleSvc(), "MdcRec/test");
587 if ( ntt ) { m_tuplet = ntt;}
588 else {
589 m_tuplet = ntupleSvc()->book ("MdcRec/test", CLID_ColumnWiseTuple, "MdcTrkRecon t particle");
590 if ( m_tuplet ) {
591 sc = m_tuplet->addItem ("layer", m_tl);
592 sc = m_tuplet->addItem ("wire", m_tw);
593 sc = m_tuplet->addItem ("phi", m_tphi);
594 sc = m_tuplet->addItem ("dphi", m_tdphi);
595 sc = m_tuplet->addItem ("dncell", m_tdncell);
596 } else {
597 log << MSG::ERROR << "Cannot book tuple MdcRec/test" << endmsg;
598 //return StatusCode::FAILURE;
599 }
600 }
601
602 //book tuple of MC truth
603 NTuplePtr ntMc(ntupleSvc(), "MdcMc/mcTk");
604 if ( ntMc ) { m_tupleMc = ntMc;}
605 else {
606 m_tupleMc = ntupleSvc()->book ("MdcMc/mcTk", CLID_ColumnWiseTuple, "MdcTrkRecon Mc particle");
607 if ( m_tupleMc ) {
608 sc = m_tupleMc->addItem ("evtNo", m_tMcEvtNo);
609 sc = m_tupleMc->addItem ("nDigi", m_tMcNTk);
610 } else {
611 log << MSG::ERROR << "Cannot book tuple MdcMc/mcTk" << endmsg;
612 //return StatusCode::FAILURE;
613 }
614 }
615
616 //book tuple of MC truth
617 NTuplePtr ntMcHit(ntupleSvc(), "MdcMc/mcHit");
618 if ( ntMcHit ) { m_tupleMcHit = ntMcHit;}
619 else {
620 m_tupleMcHit = ntupleSvc()->book ("MdcMc/mcHit", CLID_ColumnWiseTuple, "MdcTrkRecon Mc hit");
621 if ( m_tupleMcHit ) {
622 sc = m_tupleMcHit->addItem ("tkId", m_tMcTkId);
623 sc = m_tupleMcHit->addItem ("pid", m_tMcPid);
624 sc = m_tupleMcHit->addItem ("px", m_tMcPx);
625 sc = m_tupleMcHit->addItem ("py", m_tMcPy);
626 sc = m_tupleMcHit->addItem ("pz", m_tMcPz);
627 sc = m_tupleMcHit->addItem ("d0", m_tMcD0);
628 sc = m_tupleMcHit->addItem ("z0", m_tMcZ0);
629 sc = m_tupleMcHit->addItem ("theta",m_tMcTheta);
630 sc = m_tupleMcHit->addItem ("fiterm",m_tMcFiTerm);
631 sc = m_tupleMcHit->addItem ("nMcHit", m_tMcNHit, 0, 6796);
632 sc = m_tupleMcHit->addIndexedItem ("layer",m_tMcNHit, m_tMcLayer);
633 sc = m_tupleMcHit->addIndexedItem ("wire", m_tMcNHit, m_tMcWire);
634 sc = m_tupleMcHit->addIndexedItem ("drift",m_tMcNHit, m_tMcDrift);
635 sc = m_tupleMcHit->addIndexedItem ("x", m_tMcNHit, m_tMcX);
636 sc = m_tupleMcHit->addIndexedItem ("y", m_tMcNHit, m_tMcY);
637 sc = m_tupleMcHit->addIndexedItem ("z", m_tMcNHit, m_tMcZ);
638 sc = m_tupleMcHit->addIndexedItem ("lr", m_tMcNHit, m_tMcLR);
639 } else {
640 log << MSG::ERROR << "Cannot book tuple MdcMc/mcHit" << endmsg;
641 //return StatusCode::FAILURE;
642 }
643 }
644 //book tuple of recnstruction results
645 NTuplePtr nt1(ntupleSvc(), "MdcRec/rec");
646 if ( nt1 ) { m_tuple1 = nt1;}
647 else {
648 m_tuple1 = ntupleSvc()->book ("MdcRec/rec", CLID_ColumnWiseTuple, "MdcTrkRecon results");
649 if ( m_tuple1 ) {
650 sc = m_tuple1->addItem ("t0", m_t0);
651 sc = m_tuple1->addItem ("t0Stat", m_t0Stat);
652 sc = m_tuple1->addItem ("t0Truth", m_t0Truth);
653
654 sc = m_tuple1->addItem ("nTdsTk", m_nTdsTk);
655 sc = m_tuple1->addItem ("nMcTk", m_nMcTk);
656
657 sc = m_tuple1->addItem ("p", m_p);
658 sc = m_tuple1->addItem ("pt", m_pt);
659 sc = m_tuple1->addItem ("nSlay", m_nSlay);
660 sc = m_tuple1->addItem ("pz", m_pz);
661 sc = m_tuple1->addItem ("d0", m_d0);
662 sc = m_tuple1->addItem ("phi0", m_phi0);
663 sc = m_tuple1->addItem ("cpa", m_cpa);
664 sc = m_tuple1->addItem ("z0", m_z0);
665 sc = m_tuple1->addItem ("tanl", m_tanl);
666 sc = m_tuple1->addItem ("q", m_q);
667 sc = m_tuple1->addItem ("pocax", m_pocax);
668 sc = m_tuple1->addItem ("pocay", m_pocay);
669 sc = m_tuple1->addItem ("pocaz", m_pocaz);
670
671 sc = m_tuple1->addItem ("evtNo", m_evtNo);
672 sc = m_tuple1->addItem ("nSt", m_nSt);
673 sc = m_tuple1->addItem ("nAct", m_nAct);
674 sc = m_tuple1->addItem ("nDof", m_nDof);
675 sc = m_tuple1->addItem ("chi2", m_chi2);
676 sc = m_tuple1->addItem ("tkId", m_tkId);
677 sc = m_tuple1->addItem ("nHit", m_nHit, 0, 6796);
678 sc = m_tuple1->addIndexedItem ("resid", m_nHit, m_resid);
679 sc = m_tuple1->addIndexedItem ("sigma", m_nHit, m_sigma);
680 sc = m_tuple1->addIndexedItem ("driftD", m_nHit, m_driftD);
681 sc = m_tuple1->addIndexedItem ("driftT", m_nHit, m_driftT);
682 sc = m_tuple1->addIndexedItem ("doca", m_nHit, m_doca);
683 sc = m_tuple1->addIndexedItem ("entra", m_nHit, m_entra);
684 sc = m_tuple1->addIndexedItem ("ambig", m_nHit, m_ambig);
685 sc = m_tuple1->addIndexedItem ("fltLen", m_nHit, m_fltLen);
686 sc = m_tuple1->addIndexedItem ("tof", m_nHit, m_tof);
687 sc = m_tuple1->addIndexedItem ("act", m_nHit, m_act);
688 sc = m_tuple1->addIndexedItem ("tdc", m_nHit, m_tdc);
689 sc = m_tuple1->addIndexedItem ("adc", m_nHit, m_adc);
690 sc = m_tuple1->addIndexedItem ("layer", m_nHit, m_layer);
691 sc = m_tuple1->addIndexedItem ("wire", m_nHit, m_wire);
692 sc = m_tuple1->addIndexedItem ("x", m_nHit, m_x);
693 sc = m_tuple1->addIndexedItem ("y", m_nHit, m_y);
694 sc = m_tuple1->addIndexedItem ("z", m_nHit, m_z);
695 sc = m_tuple1->addIndexedItem ("dx", m_nHit, m_dx);
696 sc = m_tuple1->addIndexedItem ("dy", m_nHit, m_dy);
697 sc = m_tuple1->addIndexedItem ("dz", m_nHit, m_dz);
698 sc = m_tuple1->addIndexedItem ("dDriftD", m_nHit, m_dDriftD);
699 sc = m_tuple1->addIndexedItem ("dlr", m_nHit, m_dlr);
700 sc = m_tuple1->addIndexedItem ("cellWidth",m_nHit, m_cellWidth);//for pickHits tuning
701 } else {
702 log << MSG::ERROR << "Cannot book tuple MdcRec/rec" << endmsg;
703 //return StatusCode::FAILURE;
704 }
705 }
706 //book tuple of segment
707 NTuplePtr ntSeg(ntupleSvc(), "MdcSeg/seg");
708 if ( ntSeg ) { m_tupleSeg = ntSeg;}
709 else {
710 m_tupleSeg = ntupleSvc()->book ("MdcSeg/seg", CLID_ColumnWiseTuple, "MdcTrkRecon segment data");
711 if ( m_tupleSeg ) {
712 sc = m_tupleSeg->addItem ("evtNo", m_tsEvtNo);
713 sc = m_tupleSeg->addItem ("nSeg", m_tsNSeg);
714 sc = m_tupleSeg->addItem ("nDigi", m_tsNDigi, 0, 6796);
715 sc = m_tupleSeg->addIndexedItem ("layer", m_tsNDigi, m_tsLayer);
716 sc = m_tupleSeg->addIndexedItem ("wire", m_tsNDigi, m_tsWire);
717 sc = m_tupleSeg->addIndexedItem ("inSeg", m_tsNDigi, m_tsInSeg);
718 sc = m_tupleSeg->addIndexedItem ("mcTkId", m_tsNDigi, m_tsMcTkId);
719 } else {
720 log << MSG::ERROR << "Cannot book tuple MdcSeg/seg" << endmsg;
721 //return StatusCode::FAILURE;
722 }
723 }
724
725 //book tuple of event
726 NTuplePtr nt4(ntupleSvc(), "MdcRec/evt");
727 if ( nt4 ) { m_tupleEvt = nt4;}
728 else {
729 m_tupleEvt = ntupleSvc()->book ("MdcRec/evt", CLID_ColumnWiseTuple, "MdcTrkRecon event data");
730 if ( m_tupleEvt ) {
731 sc = m_tupleEvt->addItem ("evtNo", m_t4EvtNo);
732 sc = m_tupleEvt->addItem ("nMcTk", m_t4nMcTk );
733 sc = m_tupleEvt->addItem ("nTdsTk", m_t4nRecTk );
734 sc = m_tupleEvt->addItem ("t0", m_t4t0);
735 sc = m_tupleEvt->addItem ("t0Stat", m_t4t0Stat);
736 sc = m_tupleEvt->addItem ("t0Truth", m_t4t0Truth);
737 sc = m_tupleEvt->addItem ("nDigiUn", m_t4nDigiUnmatch);
738 sc = m_tupleEvt->addItem ("nDigi", m_t4nDigi, 0, 6796);
739 sc = m_tupleEvt->addIndexedItem ("layer", m_t4nDigi, m_t4Layer);
740 sc = m_tupleEvt->addIndexedItem ("wire", m_t4nDigi, m_t4Wire);
741 sc = m_tupleEvt->addIndexedItem ("rt", m_t4nDigi, m_t4Time);
742 sc = m_tupleEvt->addIndexedItem ("rc", m_t4nDigi, m_t4Charge);
743 sc = m_tupleEvt->addIndexedItem ("phiMid", m_t4nDigi, m_t4PhiMid);
744 sc = m_tupleEvt->addIndexedItem ("hot", m_t4nDigi, m_t4Hot);
745 //sc = m_tupleEvt->addIndexedItem ("recHit", m_t4nDigi, m_t4recHit);
746 //sc = m_tupleEvt->addIndexedItem ("rawHit", m_t4nDigi, m_t4rawHit);
747 } else {
748 log << MSG::ERROR << "Cannot book N-tuple: MdcRec/evt" << endmsg;
749 //return StatusCode::FAILURE;
750 }
751 }
752
753 //book tuple of residula for every layer
754 NTuplePtr ntCombAx(ntupleSvc(), "MdcCombAx/combAx");
755 if ( ntCombAx ) { g_tupleCombAx= ntCombAx;}
756 else {
757 g_tupleCombAx = ntupleSvc()->book ("MdcCombAx/combAx", CLID_RowWiseTuple, "MdcTrkRecon cuts in combine ax seg");
758 if ( g_tupleCombAx) {
759 sc = g_tupleCombAx->addItem ("dPhi0", g_combAxdPhi0);
760 sc = g_tupleCombAx->addItem ("dCurv", g_combAxdCurv);
761 sc = g_tupleCombAx->addItem ("sigPhi0", g_combAxSigPhi0);
762 sc = g_tupleCombAx->addItem ("sigCurv", g_combAxSigCurv);
763 sc = g_tupleCombAx->addItem ("slSeed", g_combAxSlSeed);
764 sc = g_tupleCombAx->addItem ("slTest", g_combAxSlTest);
765 sc = g_tupleCombAx->addItem ("qSeed", g_combAxQualitySeed);
766 sc = g_tupleCombAx->addItem ("qTest", g_combAxQualityTest);
767 sc = g_tupleCombAx->addItem ("pSeed", g_combAxPatSeed);
768 sc = g_tupleCombAx->addItem ("pTest", g_combAxPatTest);
769 sc = g_tupleCombAx->addItem ("nHitSeed", g_combAxNhitSeed);
770 sc = g_tupleCombAx->addItem ("nHitTest", g_combAxNhitTest);
771 sc = g_tupleCombAx->addItem ("mc", g_combAxMc);
772 sc = g_tupleCombAx->addItem ("ptTruth", g_combAxMcPt);
773 sc = g_tupleCombAx->addItem ("thetaTruth",g_combAxMcTheta);
774 sc = g_tupleCombAx->addItem ("phiTruth", g_combAxMcPhi);
775 sc = g_tupleCombAx->addItem ("ambigSeed", g_combAxMcAmbigSeed);
776 sc = g_tupleCombAx->addItem ("ambigTest", g_combAxMcAmbigTest);
777 } else {
778 log << MSG::ERROR << "Cannot book N-tuple: FILE101/combAx" << endmsg;
779 //return StatusCode::FAILURE;
780 }
781 }
782
783 //book tuple of residula for every layer
784 NTuplePtr ntFindSeg(ntupleSvc(), "MdcSeg/findSeg");
785 if ( ntFindSeg ) { g_tupleFindSeg = ntFindSeg;}
786 else {
787 g_tupleFindSeg = ntupleSvc()->book ("MdcSeg/findSeg", CLID_RowWiseTuple, "MdcTrkRecon cuts in segment finder");
788 if ( g_tupleFindSeg) {
789 sc = g_tupleFindSeg->addItem ("chi2", g_findSegChi2);
790 sc = g_tupleFindSeg->addItem ("slope", g_findSegSlope);
791 sc = g_tupleFindSeg->addItem ("intercept",g_findSegIntercept);
792 sc = g_tupleFindSeg->addItem ("chi2Refit",g_findSegChi2Refit);
793 sc = g_tupleFindSeg->addItem ("chi2Add", g_findSegChi2Add);
794 sc = g_tupleFindSeg->addItem ("pat", g_findSegPat);
795 sc = g_tupleFindSeg->addItem ("pat34", g_findSegPat34);
796 sc = g_tupleFindSeg->addItem ("nhit", g_findSegNhit);
797 sc = g_tupleFindSeg->addItem ("slayer", g_findSegSl);
798 sc = g_tupleFindSeg->addItem ("mc", g_findSegMc);
799 sc = g_tupleFindSeg->addItem ("ambig", g_findSegAmbig);
800 } else {
801 log << MSG::ERROR << "Cannot book N-tuple: FILE101/findSeg" << endmsg;
802 //return StatusCode::FAILURE;
803 }
804 }
805
806 //book tuple of event
807 NTuplePtr ntPick(ntupleSvc(), "MdcTrk/pick");
808 if ( ntPick ) { m_tuplePick = ntPick;}
809 else {
810 m_tuplePick = ntupleSvc()->book ("MdcTrk/pick", CLID_ColumnWiseTuple, "MdcTrkRecon pickHits");
811 if ( m_tuplePick ) {
812 sc = m_tuplePick->addItem ("layer", m_pickLayer);
813 sc = m_tuplePick->addItem ("nCell", m_pickNCell, 0, 288);
814 sc = m_tuplePick->addItem ("nCellWithDigi", m_pickNCellWithDigi, 0, 288);
815 sc = m_tuplePick->addIndexedItem ("wire", m_pickNCellWithDigi, m_pickWire);
816 sc = m_tuplePick->addIndexedItem ("predDrift", m_pickNCellWithDigi, m_pickPredDrift);
817 sc = m_tuplePick->addIndexedItem ("drift", m_pickNCellWithDigi, m_pickDrift);
818 sc = m_tuplePick->addIndexedItem ("driftTruth", m_pickNCellWithDigi, m_pickDriftTruth);
819 sc = m_tuplePick->addIndexedItem ("phiZOk", m_pickNCellWithDigi, m_pickPhizOk);
820 sc = m_tuplePick->addIndexedItem ("z", m_pickNCellWithDigi, m_pickZ);
821 sc = m_tuplePick->addIndexedItem ("resid", m_pickNCellWithDigi, m_pickResid);
822 sc = m_tuplePick->addIndexedItem ("sigma", m_pickNCellWithDigi, m_pickSigma);
823 sc = m_tuplePick->addIndexedItem ("phiLowCut", m_pickNCellWithDigi, m_pickPhiLowCut);
824 sc = m_tuplePick->addIndexedItem ("phiHighCut", m_pickNCellWithDigi, m_pickPhiHighCut);
825 sc = m_tuplePick->addIndexedItem ("goodDriftCut", m_pickNCellWithDigi, m_pickDriftCut);
826 sc = m_tuplePick->addIndexedItem ("mcTk", m_pickNCellWithDigi, m_pickMcTk);
827 sc = m_tuplePick->addIndexedItem ("is2d", m_pickNCellWithDigi, m_pickIs2D);
828 sc = m_tuplePick->addIndexedItem ("pt", m_pickNCellWithDigi, m_pickPt);
829 sc = m_tuplePick->addIndexedItem ("curv", m_pickNCellWithDigi, m_pickCurv);
830 } else {
831 log << MSG::ERROR << "Cannot book N-tuple: MdcTrk/pick" << endmsg;
832 //return StatusCode::FAILURE;
833 }
834 }
835
836#ifdef MDCPATREC_TIMETEST
837 //book tuple of time test
838 NTuplePtr nt5(ntupleSvc(), "MdcTime/ti");
839 if ( nt5 ) { m_tupleTime = nt5;}
840 else {
841 m_tupleTime = ntupleSvc()->book ("MdcTime/ti", CLID_RowWiseTuple, "MdcTrkRecon time");
842 if ( m_tupleTime ) {
843 sc = m_tupleTime->addItem ("eventtime", m_eventTime);
844 sc = m_tupleTime->addItem ("recTkNum", m_t5recTkNum);
845 sc = m_tupleTime->addItem ("mcTkNum", m_mcTkNum);
846 sc = m_tupleTime->addItem ("evtNo", m_t5EvtNo);
847 sc = m_tupleTime->addItem ("nHit", m_t5nHit);
848 sc = m_tupleTime->addItem ("nDigi", m_t5nDigi);
849 } else {
850 log << MSG::ERROR << "Cannot book N-tuple: FILE101/time" << endmsg;
851 //return StatusCode::FAILURE;
852 }
853 }
854 sc = service( "BesTimerSvc", m_timersvc);
855 if( sc.isFailure() ) {
856 log << MSG::WARNING << " Unable to locate BesTimerSvc" << endreq;
857 //return StatusCode::FAILURE;
858 }
859 m_timer[1] = m_timersvc->addItem("Execution");
860 m_timer[1]->propName("nExecution");
861#endif
862
863#ifdef MDCPATREC_RESLAYER
864 //book tuple of residula for every layer
865 NTuplePtr nt6(ntupleSvc(), "MdcRes/res");
866 if ( nt6 ) { g_tupleres = nt6;}
867 else {
868 g_tupleres= ntupleSvc()->book ("MdcRes/res", CLID_RowWiseTuple, "MdcTrkRecon res");
869 if ( g_tupleres ) {
870 sc = g_tupleres->addItem ("resid", g_resLayer);
871 sc = g_tupleres->addItem ("layer", g_t6Layer);
872 } else {
873 log << MSG::ERROR << "Cannot book N-tuple: FILE101/res" << endmsg;
874 return StatusCode::FAILURE;
875 }
876 }
877#endif
878
879 return sc;
880}// end of bookNTuple()
881
882void MdcTrkRecon::fillMcTruth(){
883 MsgStream log(msgSvc(), name());
884 StatusCode sc;
885 t_mcTkNum = 0;
886 for(int i=0;i<100;i++){
887 isPrimaryOfMcTk[i]=-9999;
888 pdgOfMcTk[i]=-9999;
889 }
890 double ptOfMcTk[100]={0.};
891 double thetaOfMcTk[100]={0.};
892 double phiOfMcTk[100]={0.};
893 double nDigiOfMcTk[100]={0.};
894 if (m_mcHist){
895 //------------------Retrieve MC truth MdcMcHit------------
896 SmartDataPtr<Event::MdcMcHitCol> mcMdcMcHitCol(eventSvc(),"/Event/MC/MdcMcHitCol");
897 if (!mcMdcMcHitCol) {
898 log << MSG::INFO << "Could not find MdcMcHit" << endreq;
899 }else{
900 Event::MdcMcHitCol::iterator iter_mchit = mcMdcMcHitCol->begin();
901 for (;iter_mchit != mcMdcMcHitCol->end(); iter_mchit++ ) {
902 const Identifier id= (*iter_mchit)->identify();
903 int layer = MdcID::layer(id);
904 int wire = MdcID::wire(id);
905 int iMcTk = (*iter_mchit)->getTrackIndex();
906 if(iMcTk<100&&iMcTk>0) nDigiOfMcTk[iMcTk]++;
907 mcDrift[layer][wire]= (*iter_mchit)->getDriftDistance()/10.; //drift in MC.
908 //std::cout<<" "<<__FILE__<<" mcDrift "<<mcDrift[layer][wire]<<" "<<std::endl;
909 mcX[layer][wire] = (*iter_mchit)->getPositionX()/10.;
910 mcY[layer][wire] = (*iter_mchit)->getPositionY()/10.;
911 mcZ[layer][wire] = (*iter_mchit)->getPositionZ()/10.;
912 mcLR[layer][wire] = (*iter_mchit)->getPositionFlag();
913 if (mcLR[layer][wire] == 0) mcLR[layer][wire] = -1;
914 t_nHitInTk[iMcTk]++;
915 haveDigiAmbig[layer][wire] = mcLR[layer][wire];
916 haveDigiDrift[layer][wire] = mcDrift[layer][wire];
917 }
918 }
919
920 //------------------get event start time truth-----------
921 t_t0Truth = -10.;
922 SmartDataPtr<Event::McParticleCol> mcParticleCol(eventSvc(),"/Event/MC/McParticleCol");
923 if(!mcParticleCol){
924 log << MSG::INFO << "Could not retrieve McParticelCol" << endreq;
925 }else {
926 int nMcTk = 0;
927 Event::McParticleCol::iterator it = mcParticleCol->begin();
928 for (;it != mcParticleCol->end(); it++){
929 //int pdg_code = (*it)->particleProperty();
930 //if ((fabs(pdg_code)!=211) || (fabs(pdg_code)!=13)) continue;
931 t_mcTkNum++;
932 nMcTk++;
933 }
934 it = mcParticleCol->begin();
935 for (;it != mcParticleCol->end(); it++){
936 int tkId = (*it)->trackIndex();
937 if ((*it)->primaryParticle()){
938 t_t0Truth = (*it)->initialPosition().t();
939 isPrimaryOfMcTk[tkId] = 1;
940 }else{
941 isPrimaryOfMcTk[tkId] = 0;
942 //continue;
943 }
944 //fill charged particle
945 int pdg_code = (*it)->particleProperty();
946 pdgOfMcTk[tkId] = pdg_code;
947 //std::cout<<" tkId "<<tkId<<" pdg_code "<<pdg_code<<" "<<std::endl;
948 //FIXME skip charged track and track outside MDC
949 //if ((fabs(pdg_code)!=211) || (fabs(pdg_code)!=13)) continue;
950 const CLHEP::Hep3Vector& true_mom = (*it)->initialFourMomentum().vect();
951 const CLHEP::HepLorentzVector& posIni = (*it)->initialPosition();
952 const CLHEP::HepLorentzVector& posFin = (*it)->finalPosition();
953 if(tkId<100&&tkId>=0) {
954 ptOfMcTk[tkId] = true_mom.perp();
955 thetaOfMcTk[tkId] = (*it)->initialFourMomentum().theta();
956 phiOfMcTk[tkId] = (*it)->initialFourMomentum().phi();
957 }
958
959 if ( m_tupleMcHit){
960 m_tMcPx = true_mom.x();
961 m_tMcPy = true_mom.y();
962 m_tMcPz = true_mom.z();
963 m_tMcD0 = posIni.perp()/10.;
964 m_tMcZ0 = posIni.z()/10.;
965 m_tMcTheta = (*it)->initialFourMomentum().theta();
966 m_tMcFiTerm= posFin.phi();
967 m_tMcPid = pdg_code;
968 m_tMcTkId = tkId;
969 m_tMcNHit = t_nHitInTk[tkId];
970 m_tupleMcHit->write();
971 }
972 }//end of loop mcParticleCol
973 if(m_tupleMc) {
974 m_tMcNTk= nMcTk;
975 m_tMcEvtNo = t_eventNo;
976 m_tupleMc->write();
977 }
978 }
979 }//end of m_mcHist
980
981 uint32_t getDigiFlag = 0;
982 getDigiFlag += m_maxMdcDigi;
983 if(m_dropHot) getDigiFlag |= MdcRawDataProvider::b_dropHot;
984 if(m_keepBadTdc) getDigiFlag |= MdcRawDataProvider::b_keepBadTdc;
985 if(m_keepUnmatch) getDigiFlag |= MdcRawDataProvider::b_keepUnmatch;
986 MdcDigiVec mdcDigiVec = m_rawDataProviderSvc->getMdcDigiVec(getDigiFlag);
987 if ((int)mdcDigiVec.size()<m_minMdcDigi){
988 std::cout << " Skip this event for MdcDigiVec.size() < "<<m_minMdcDigi << endl;
989 }
990
991 if (mdcDigiVec.size()<=0) return;
992 //fill raw digi and t0 status
993 MdcDigiCol::iterator iter = mdcDigiVec.begin();
994 for (;iter!= mdcDigiVec.end(); iter++ ) {
995 long l = MdcID::layer((*iter)->identify());
996 long w = MdcID::wire((*iter)->identify());
997 int tkId = (*iter)->getTrackIndex();
998 haveDigiTk[l][w]= tkId;
999 //std::cout<<" l "<<l<<" w "<<w<<" tk "<<tkId<<std::endl;
1000 haveDigiPt[l][w] = ptOfMcTk[(*iter)->getTrackIndex()];
1001 haveDigiTheta[l][w] = thetaOfMcTk[(*iter)->getTrackIndex()];
1002 haveDigiPhi[l][w] = phiOfMcTk[(*iter)->getTrackIndex()];
1003 if(m_tupleWireEff){
1004 m_we_tkId = tkId;
1005 if(tkId>=0){
1006 if(tkId>=1000) tkId-=1000;
1007 m_we_pdg = pdgOfMcTk[tkId];
1008 m_we_primary = isPrimaryOfMcTk[tkId];
1009 //if(pdgOfMcTk[tkId]==-22) cout<<" primaryParticle ? "<< isPrimaryOfMcTk[tkId]<< " tkId "<<tkId <<" layer "<<l<<" wire "<<w<<" pdg="<<pdgOfMcTk[tkId]<<endl;
1010 }else{
1011 m_we_pdg = -999;
1012 m_we_primary = -999;
1013 }
1014 m_we_layer = l;
1015 m_we_wire = w;
1016 int gwire = w+Constants::nWireBeforeLayer[l];
1017 m_we_gwire = gwire;
1018 m_we_poisoned = hitPoisoned[l][w];
1019 if(hitInSegList[l][w]>0) m_we_seg = 1;
1020 else m_we_seg = 0;
1021 if(recHitMap[l][w]>0) m_we_track = 1;
1022 else m_we_track = 0;
1023 if(l==35&&tkId>=0&&abs(pdgOfMcTk[tkId])==11&&hitInSegList[l][w]<=0) {
1024 cout<<"layer "<<l<<" cell "<<w<<" gwire "<<gwire<<" inseg "<<hitInSegList[l][w]<<endl;
1025 }
1026 m_we_pt = ptOfMcTk[tkId];
1027 m_we_theta = thetaOfMcTk[tkId];
1028 m_we_phi = phiOfMcTk[tkId];
1029 m_tupleWireEff->write();
1030 }
1031 }
1032}//end of fillMcTruth()
1033
1034void MdcTrkRecon::fillSegList() {
1035 if (2 == m_flags.segPar.lPrint) {
1036 std::cout << "print after Segment finding " << std::endl;
1037 }
1038 if (!m_flags.lHist) return;
1039 // Fill hits of every layer after segment finding
1040 for (int ii=0;ii<43;ii++){
1041 for (int jj=0;jj<288;jj++){ hitInSegList[ii][jj] = 0; }
1042 }
1043 int nSeg=0;
1044 for (int i = 0; i < _gm->nSuper(); i++) {
1045 MdcSeg* aSeg = (MdcSeg *) m_segs->oneList(i)->first();
1046 while (aSeg != NULL) {
1047 nSeg++;
1048 for (int ihit=0 ; ihit< aSeg->nHit() ; ihit++){
1049 const MdcHit* hit = aSeg->hit(ihit)->mdcHit();
1050 int layer = hit->layernumber();
1051 int wire = hit->wirenumber();
1052 hitInSegList[layer][wire]++;
1053 }
1054 aSeg = (MdcSeg *) aSeg->next();
1055 }
1056 }//end for slayer
1057 if (!m_tupleSeg) return;
1058 m_tsEvtNo = t_eventNo;
1059 m_tsNDigi = t_nDigi;
1060 int iDigi=0;
1061 for (int ii=0;ii<43;ii++){
1062 for (int jj=0;jj<288;jj++){
1063 if (haveDigiTk[ii][jj] > -2){
1064 m_tsLayer[iDigi] = ii;
1065 m_tsWire[iDigi] = jj;
1066 m_tsInSeg[iDigi] = hitInSegList[ii][jj];
1067 m_tsMcTkId[iDigi] = haveDigiTk[ii][jj];
1068 iDigi++;
1069 }
1070 }
1071 }
1072 m_tsNSeg=nSeg;
1073 m_tupleSeg->write();
1074}//end of fillSegList
1075
1076void MdcTrkRecon::dumpTdsTrack(RecMdcTrackCol* trackList){
1077 std::cout<< "tksize = "<<trackList->size() << std::endl;//yzhang debug
1078 RecMdcTrackCol::iterator it = trackList->begin();
1079 for (;it!= trackList->end();it++){
1080 RecMdcTrack *tk = *it;
1081 std::cout<< endl<<"//====RecMdcTrack "<<tk->trackId()<<"====:" << std::endl;
1082 cout <<" d0 "<<tk->helix(0)
1083 <<" phi0 "<<tk->helix(1)
1084 <<" cpa "<<tk->helix(2)
1085 <<" z0 "<<tk->helix(3)
1086 <<" tanl "<<tk->helix(4)
1087 <<endl;
1088 std::cout<<" q "<<tk->charge()
1089 <<" theta "<<tk->theta()
1090 <<" phi "<<tk->phi()
1091 <<" x0 "<<tk->x()
1092 <<" y0 "<<tk->y()
1093 <<" z0 "<<tk->z()
1094 <<" r0 "<<tk->r()
1095 <<endl;
1096 std::cout <<" p "<<tk->p()
1097 <<" pt "<<tk->pxy()
1098 <<" px "<<tk->px()
1099 <<" py "<<tk->py()
1100 <<" pz "<<tk->pz()
1101 <<endl;
1102 std::cout<<" tkStat "<<tk->stat()
1103 <<" chi2 "<<tk->chi2()
1104 <<" ndof "<<tk->ndof()
1105 <<" nhit "<<tk->getNhits()
1106 <<" nst "<<tk->nster()
1107 <<endl;
1108 std::cout<< "errmat mat " << std::endl;
1109 std::cout<< tk->err()<<std::endl;
1110 //std::cout<< "errmat array " << std::endl;
1111 //for (int i=0; i<15; i++){ std::cout<< " "<<tk->err(i); }
1112 //std::cout<< " " << std::endl;
1113
1114 int nhits = tk->getVecHits().size();
1115 std::cout<<nhits <<" Hits: " << std::endl;
1116 for(int ii=0; ii <nhits ; ii++){
1117 Identifier id(tk->getVecHits()[ii]->getMdcId());
1118 int layer = MdcID::layer(id);
1119 int wire = MdcID::wire(id);
1120 cout<<"("<< layer <<","<<wire<<","<<tk->getVecHits()[ii]->getStat()
1121 <<",lr:"<<tk->getVecHits()[ii]->getFlagLR()<<") ";
1122 }//end of hit list
1123 std::cout << " "<< std::endl;
1124 }//end of tk list
1125 std::cout << " "<< std::endl;
1126}//end of dumpTdsTrack
1127
1128void MdcTrkRecon::fillTrackList(){
1129 if (m_flags.debugFlag()>1) {
1130 std::cout << "=======Print after Track finding=======" << std::endl;
1131 m_tracks->plot();
1132 }
1133
1134 if (!m_flags.lHist) return;
1135 //----------- fill hit -----------
1136 int nTk = (*m_tracks).nTrack();
1137 for (int itrack = 0; itrack < nTk; itrack++) {
1138 MdcTrack *atrack = (*m_tracks)[itrack];
1139 if (atrack==NULL) continue;
1140 const TrkFit* fitResult = atrack->track().fitResult();
1141 if (fitResult == 0) continue;//check the fit worked
1142
1143 //for fill ntuples
1144 int nSt=0; //int nSeg=0;
1145 int seg[11] = {0}; int segme;
1146 //-----------hit list-------------
1147 const TrkHitList* hitlist = atrack->track().hits();
1148 TrkHitList::hot_iterator hot = hitlist->begin();
1149 int layerCount[43]={0};
1150 for(int iii=0;iii<43;iii++){layerCount[iii]=0;}
1151 int i=0;
1152 for (;hot!= hitlist->end();hot++){
1153 const MdcRecoHitOnTrack* recoHot;
1154 recoHot = dynamic_cast<const MdcRecoHitOnTrack*> (&(*hot));
1155 int layer = recoHot->mdcHit()->layernumber();
1156 int wire = recoHot->mdcHit()->wirenumber();
1157 layerCount[layer]++;
1158 recHitMap[layer][wire]++;
1159 //for n seg
1160 segme=0;
1161 if ( layer >0 ) {segme=(layer-1)/4;}
1162 seg[segme]++;
1163 if (recoHot->layer()->view()) { ++nSt; }
1164
1165 if(!m_tuple1) continue;
1166 m_layer[i] = layer;
1167 m_wire[i] = wire;
1168 m_ambig[i] = recoHot->wireAmbig();// wire ambig
1169 //fltLen
1170 double fltLen = recoHot->fltLen();
1171 m_fltLen[i] = fltLen;
1172 double tof = recoHot->getParentRep()->arrivalTime(fltLen);
1173 //position
1174 HepPoint3D pos; Hep3Vector dir;
1175 recoHot->getParentRep()->traj().getInfo(fltLen,pos,dir);
1176 m_x[i] = pos.x();
1177 m_y[i] = pos.y();
1178 m_z[i] = pos.z();
1179 m_driftT[i] = recoHot->mdcHit()->driftTime(tof,pos.z());
1180 m_tof[i] = tof;
1181 m_driftD[i] = recoHot->drift();
1182 m_sigma[i] = recoHot->hitRms();
1183 m_tdc[i] = recoHot->rawTime();
1184 m_adc[i] = recoHot->mdcHit()->charge();
1185 m_doca[i] = recoHot->dcaToWire();//sign w.r.t. dirft() FIXME
1186 m_entra[i] = recoHot->entranceAngle();
1187 m_act[i] = recoHot->isActive();
1188 //diff with truth
1189 m_dx[i] = m_x[i] - mcX[layer][wire];
1190 m_dy[i] = m_y[i] - mcY[layer][wire];
1191 m_dz[i] = m_z[i] - mcZ[layer][wire];
1192 m_dDriftD[i] = m_driftD[i] - mcDrift[layer][wire];
1193 m_dlr[i] = m_ambig[i] - mcLR[layer][wire];
1194 //yzhang for pickHits debug
1195 m_cellWidth[i] = Constants::twoPi*_gm->Layer(layer)->rMid()/_gm->Layer(layer)->nWires();
1196 //zhangy
1197 //resid
1198 double res=999.,rese=999.;
1199 if (recoHot->resid(res,rese,false)){
1200 }else{}
1201 m_resid[i] = res;
1202 i++;
1203 }// end fill hit
1204 int nSlay=0;
1205 for(int i=0;i<11;i++){
1206 if (seg[i]>0) nSlay++;
1207 }
1208 if(m_tuple1){
1209 //------------fill track result-------------
1210 m_tkId = itrack;
1211 //track parameters at closest approach to beamline
1212 double fltLenPoca = 0.0;
1213 TrkExchangePar helix = fitResult->helix(fltLenPoca);
1214 m_q = fitResult->charge();
1215 m_phi0 = helix.phi0();
1216 m_tanl = helix.tanDip();
1217 m_z0 = helix.z0();
1218 m_d0 = helix.d0();
1219 m_pt = fitResult->pt();
1220 m_nSlay = nSlay;
1221 if (m_pt > 0.00001){
1222 m_cpa = fitResult->charge()/fitResult->pt();
1223 }
1224 //momenta and position
1225 CLHEP::Hep3Vector p1 = fitResult->momentum(fltLenPoca);
1226 double px,py,pz,pxy;
1227 pxy = fitResult->pt();
1228 px = p1.x();
1229 py = p1.y();
1230 pz = p1.z();
1231 m_p = p1.mag();
1232 m_pz = pz;
1233 HepPoint3D poca = fitResult->position(fltLenPoca);
1234 m_pocax = poca.x();
1235 m_pocay = poca.y();
1236 m_pocaz = poca.z();
1237 m_nAct = fitResult->nActive();
1238 m_nHit = hitlist->nHit();
1239 m_nDof = fitResult->nDof();
1240 m_nSt = nSt;
1241 m_chi2 = fitResult->chisq();
1242 //for (int l=0;l<43;l++) m_layerCount[l] = layerCount[l];
1243 m_t0 = t_t0;
1244 m_t0Stat = t_t0Stat;
1245 m_t0Truth = t_t0Truth;
1246 m_nTdsTk = t_nRecTk;
1247 m_nMcTk = t_mcTkNum;
1248 m_evtNo = t_eventNo;
1249 m_tuple1->write();
1250 }
1251 }//end of loop rec trk list
1252}//end of fillTrackList
1253
1254
1255void MdcTrkRecon::fillEvent(){
1256 if (m_tupleEvt!=NULL){
1257
1258 uint32_t getDigiFlag = 0;
1259 getDigiFlag += m_maxMdcDigi;
1260 if(m_dropHot) getDigiFlag |= MdcRawDataProvider::b_dropHot;
1261 if(m_keepBadTdc) getDigiFlag |= MdcRawDataProvider::b_keepBadTdc;
1262 if(m_keepUnmatch) getDigiFlag |= MdcRawDataProvider::b_keepUnmatch;
1263 MdcDigiVec mdcDigiVec = m_rawDataProviderSvc->getMdcDigiVec(getDigiFlag);
1264 if ((int)mdcDigiVec.size()<m_minMdcDigi){
1265 std::cout << " Skip this event for MdcDigiVec.size() < "<<m_minMdcDigi << endl;
1266 }
1267
1268 m_t4nDigi = mdcDigiVec.size();
1269 int iDigi=0;
1270 MdcDigiCol::iterator iter = mdcDigiVec.begin();
1271 for (;iDigi<m_t4nDigi; iter++ ) {
1272 double tdc = RawDataUtil::MdcTime((*iter)->getTimeChannel());
1273 double adc = RawDataUtil::MdcCharge((*iter)->getChargeChannel());
1274 m_t4Time[iDigi] = tdc;
1275 m_t4Charge[iDigi] = adc;
1276 long l = MdcID::layer((*iter)->identify());
1277 long w = MdcID::wire((*iter)->identify());
1278 m_t4Layer[iDigi] = l;
1279 m_t4Wire[iDigi] = w;
1280 m_t4PhiMid[iDigi] = _gm->Layer(l)->phiWire(w);
1281 m_t4Hot[iDigi] = recHitMap[l][w];
1282 iDigi++;
1283 }
1284 m_t4t0 = t_t0;
1285 m_t4t0Stat = t_t0Stat;
1286 m_t4t0Truth = t_t0Truth;
1287 m_t4EvtNo = t_eventNo;
1288 m_t4nRecTk = t_nRecTk;
1289 SmartDataPtr<MdcDigiCol> mdcDigiCol(eventSvc(),"/Event/Digi/MdcDigiCol");
1290 if (mdcDigiCol){
1291 m_t4nDigiUnmatch = mdcDigiCol->size();
1292 }
1293 m_t4nMcTk= t_mcTkNum;
1294 m_tupleEvt->write();
1295 }
1296}//end of fillEvent
1297
1298void MdcTrkRecon::dumpDigi(){
1299 uint32_t getDigiFlag = 0;
1300 getDigiFlag += m_maxMdcDigi;
1301 if(m_dropHot) getDigiFlag |= MdcRawDataProvider::b_dropHot;
1302 if(m_keepBadTdc) getDigiFlag |= MdcRawDataProvider::b_keepBadTdc;
1303 if(m_keepUnmatch) getDigiFlag |= MdcRawDataProvider::b_keepUnmatch;
1304 MdcDigiVec mdcDigiVec = m_rawDataProviderSvc->getMdcDigiVec(getDigiFlag);
1305 MdcDigiCol::iterator iter = mdcDigiVec.begin();
1306 std::cout<<name()<<" dump MdcDigiVec, nDigi="<<mdcDigiVec.size()<<std::endl;
1307 for (int iDigi=0;iter!= mdcDigiVec.end(); iter++,iDigi++ ) {
1308 long l = MdcID::layer((*iter)->identify());
1309 long w = MdcID::wire((*iter)->identify());
1310 int tkTruth = (*iter)->getTrackIndex();
1311 double tdc = RawDataUtil::MdcTime((*iter)->getTimeChannel());
1312 double adc = RawDataUtil::MdcCharge((*iter)->getChargeChannel());
1313 std::cout<<"("<<l<<","<<w<<";"<<tkTruth<<",t "<<tdc<<",c "<<adc<<")";
1314 if(iDigi%4==0) std::cout<<std::endl;
1315 }
1316 std::cout<<std::endl;
1317}
w
double w
Definition: EvtPhokharaDef.hh:26
iter
EvtStreamInputIterator< typename Generator::result_type > iter(Generator gen, int N=0)
Definition: EvtStreamInputIterator.hh:100
RecMdcHitCol
ObjectVector< RecMdcHit > RecMdcHitCol
Definition: InstallArea/include/MdcRecEvent/MdcRecEvent/RecMdcHit.h:99
RecMdcTrackCol
ObjectVector< RecMdcTrack > RecMdcTrackCol
Definition: InstallArea/include/MdcRecEvent/MdcRecEvent/RecMdcTrack.h:79
g_combAxQualitySeed
NTuple::Item< double > g_combAxQualitySeed
Definition: InstallArea/include/MdcTrkRecon/MdcTrkRecon/MdcHistItem.h:156
g_combAxMcAmbigTest
NTuple::Item< double > g_combAxMcAmbigTest
Definition: InstallArea/include/MdcTrkRecon/MdcTrkRecon/MdcHistItem.h:167
m_pickPhiLowCut
NTuple::Array< double > m_pickPhiLowCut
Definition: InstallArea/include/MdcTrkRecon/MdcTrkRecon/MdcHistItem.h:224
g_findSegChi2Refit
NTuple::Item< double > g_findSegChi2Refit
Definition: InstallArea/include/MdcTrkRecon/MdcTrkRecon/MdcHistItem.h:201
m_pickDriftCut
NTuple::Array< double > m_pickDriftCut
Definition: InstallArea/include/MdcTrkRecon/MdcTrkRecon/MdcHistItem.h:226
g_findSegIntercept
NTuple::Item< double > g_findSegIntercept
Definition: InstallArea/include/MdcTrkRecon/MdcTrkRecon/MdcHistItem.h:199
g_combAxQualityTest
NTuple::Item< double > g_combAxQualityTest
Definition: InstallArea/include/MdcTrkRecon/MdcTrkRecon/MdcHistItem.h:157
m_pickDriftTruth
NTuple::Array< double > m_pickDriftTruth
Definition: InstallArea/include/MdcTrkRecon/MdcTrkRecon/MdcHistItem.h:217
m_pickPredDrift
NTuple::Array< double > m_pickPredDrift
Definition: InstallArea/include/MdcTrkRecon/MdcTrkRecon/MdcHistItem.h:215
g_combAxMcAmbigSeed
NTuple::Item< double > g_combAxMcAmbigSeed
Definition: InstallArea/include/MdcTrkRecon/MdcTrkRecon/MdcHistItem.h:166
m_pickPhiHighCut
NTuple::Array< double > m_pickPhiHighCut
Definition: InstallArea/include/MdcTrkRecon/MdcTrkRecon/MdcHistItem.h:225
DstMdcTrack::err
const HepSymMatrix err() const
Definition: DstMdcTrack.cxx:155
DstMdcTrack::helix
const HepVector helix() const
......
Definition: DstMdcTrack.cxx:151
MdcDetector::Layer
const MdcLayer * Layer(unsigned id) const
Definition: InstallArea/include/MdcGeom/MdcGeom/MdcDetector.h:33
MdcDetector::instance
static MdcDetector * instance()
Definition: MdcDetector.cxx:21
MdcFlagHold::setDebug
void setDebug(int debugFlag)
Definition: MdcFlagHold.cxx:18
MdcFlagHold::readPar
void readPar(std::string inname)
Definition: MdcFlagHold.cxx:35
MdcHitOnTrack::wireAmbig
int wireAmbig() const
Definition: MdcHitOnTrack.cxx:415
MdcHitOnTrack::rawTime
double rawTime() const
Definition: MdcHitOnTrack.cxx:372
MdcHitOnTrack::dcaToWire
double dcaToWire() const
Definition: MdcHitOnTrack.cxx:123
MdcHitOnTrack::entranceAngle
double entranceAngle() const
Definition: MdcHitOnTrack.cxx:172
MdcHitOnTrack::layer
const MdcLayer * layer() const
Definition: MdcHitOnTrack.cxx:354
MdcHitUse::mdcHit
const MdcHit * mdcHit() const
Definition: MdcHitUse.cxx:69
MdcHit::driftTime
double driftTime(double tof, double z) const
Definition: MdcHit.cxx:142
MdcID::layer
static int layer(const Identifier &id)
Values of different levels (failure returns 0)
Definition: MdcID.cxx:49
MdcID::wire
static int wire(const Identifier &id)
Definition: MdcID.cxx:54
MdcLayer::phiWire
double phiWire(int cell) const
Definition: MdcLayer.cxx:87
MdcPrintSvc::printRecMdcTrackCol
void printRecMdcTrackCol() const
Definition: MdcPrintSvc.cxx:109
MdcRecoHitOnTrack::mdcHit
const MdcHit * mdcHit() const
Definition: MdcRecoHitOnTrack.cxx:38
MdcSeg::nHit
int nHit() const
Definition: MdcSeg.cxx:368
MdcTrkRecon::bookNTuple
StatusCode bookNTuple()
Definition: MdcTrkRecon.cxx:527
MdcTrkRecon::finalize
StatusCode finalize()
Definition: MdcTrkRecon.cxx:516
MdcTrkRecon::dumpTdsTrack
void dumpTdsTrack(RecMdcTrackCol *trackList)
Definition: MdcTrkRecon.cxx:1076
MdcTrkRecon::fillEvent
void fillEvent()
Definition: MdcTrkRecon.cxx:1255
MdcTrkRecon::execute
StatusCode execute()
Definition: MdcTrkRecon.cxx:217
MdcTrkRecon::fillTrackList
void fillTrackList()
Definition: MdcTrkRecon.cxx:1128
MdcTrkRecon::initialize
StatusCode initialize()
Definition: MdcTrkRecon.cxx:139
MdcTrkRecon::beginRun
StatusCode beginRun()
Definition: MdcTrkRecon.cxx:123
MdcTrkRecon::MdcTrkRecon
MdcTrkRecon(const std::string &name, ISvcLocator *pSvcLocator)
Definition: MdcTrkRecon.cxx:71
MdcTrkRecon::fillSegList
void fillSegList()
Definition: MdcTrkRecon.cxx:1034
MdcTrkRecon::fillMcTruth
void fillMcTruth()
Definition: MdcTrkRecon.cxx:882
Pdt::readMCppTable
static void readMCppTable(std::string filenm)
RawDataUtil::MdcCharge
static double MdcCharge(int chargeChannel)
Definition: Event/RawEvent/RawEvent-00-03-19/RawEvent/RawDataUtil.h:11
Trajectory::getInfo
virtual void getInfo(double fltLen, HepPoint3D &pos, Hep3Vector &direction) const =0
TrkAbsFit::momentum
virtual Hep3Vector momentum(double fltL=0.) const =0
TrkAbsFit::pt
virtual double pt(double fltL=0.) const =0
TrkAbsFit::chisq
virtual double chisq() const =0
TrkAbsFit::charge
virtual int charge() const =0
TrkAbsFit::nDof
virtual int nDof() const =0
TrkAbsFit::traj
virtual const TrkDifTraj & traj() const =0
TrkAbsFit::position
virtual HepPoint3D position(double fltL) const =0
TrkFit::nActive
virtual int nActive() const =0
TrkFit::helix
virtual TrkExchangePar helix(double fltL) const =0
TrkHitOnTrk::resid
double resid(bool exclude=false) const
Definition: TrkHitOnTrk.cxx:191
TrkRecoTrk::fitResult
const TrkFit * fitResult() const
Definition: TrkRecoTrk.cxx:387
TrkRep::arrivalTime
virtual double arrivalTime(double fltL) const
Definition: TrkRep.cxx:192
nhits
int nhits
Definition: prepare_data.cxx:29