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