BOSS 7.0.7
BESIII Offline Software System
Loading...
Searching...
No Matches
MdcxTrackFinder.cxx
Go to the documentation of this file.
1//--------------------------------------------------------------------------
2// File and Version Information:
3// $Id: MdcxTrackFinder.cxx,v 1.50 2017/08/18 06:51:52 zhangy Exp $
4//
5// Description:
6// Class MdcxTrackFinder
7//
8// Environment:
9// Software developed for the BaBar Detector at the SLAC B-Factory.
10//
11// Author List:
12// Steve Wagner Original Author
13// Zhang Yao([email protected])
14//
15// Copyright Information:
16// Copyright (C) 1997 BEPCII
17//
18// History:
19// Migration for BESIII MDC
20//
21//------------------------------------------------------------------------
22
23//-----------------------
24// This Class's Header --
25//-----------------------
27
28//-------------
29// C Headers --
30//-------------
31#include <assert.h>
32#include <time.h>
33#include <math.h>
34
35//---------------
36// C++ Headers --
37//---------------
38#include <iostream>
39
40//-------------------------------
41// Collaborating Class Headers --
42//-------------------------------
43#include "GaudiKernel/MsgStream.h"
44#include "GaudiKernel/AlgFactory.h"
45#include "BField/BField.h"
46#include "MdcRawEvent/MdcDigi.h"
47#include "MdcGeom/MdcDetector.h"
49
50#include "MdcxReco/MdcxHit.h"
51#include "MdcxReco/MdcxHits.h"
53#include "MdcxReco/MdcxSeg.h"
54#include "MdcxReco/MdcxHel.h"
57#include "MdcxReco/Mdcxprobab.h"
58#include "CLHEP/Alist/AIterator.h"
60#include "TrkBase/TrkErrCode.h"
61#include "TrkBase/TrkRecoTrk.h"
62#include "TrkBase/TrkFit.h"
63#include "TrkBase/TrkHitList.h"
72#include "Identifier/MdcID.h"
75#include "McTruth/McParticle.h"
76#include "McTruth/MdcMcHit.h"
81#include "MdcRecoUtil/Pdt.h"
83//debug
84#include "TrkBase/TrkHotList.h"
85#include "GaudiKernel/INTupleSvc.h"
86#include "GaudiKernel/NTuple.h"
88#include "MdcGeom/Constants.h"
91#include "TrigEvent/TrigData.h"
92
96
97
98
99using std::cout;
100using std::endl;
101
102extern double MdcTrkReconCut_helix_fit[43];
103//----------------
104// Constructors --
105//----------------
106
107MdcxTrackFinder::MdcxTrackFinder(const std::string& name, ISvcLocator* pSvcLocator) :
108 Algorithm(name, pSvcLocator), m_mdcCalibFunSvc(0)
109{
110 //input
111 declareProperty("pdtFile", m_pdtFile = "pdt.table");
112 //debug control
113 declareProperty("debug", m_debug= 0);
114 declareProperty("hist", m_hist= 0);
115 declareProperty("mcHist", m_mcHist = false);
116 //cuts and control
117 declareProperty("cresol", m_cresol = 0.013);
118
119 declareProperty("getDigiFlag", m_getDigiFlag = 0);
120 declareProperty("maxMdcDigi", m_maxMdcDigi= 0);
121 declareProperty("keepBadTdc", m_keepBadTdc= 0);
122 declareProperty("dropHot", m_dropHot= 0);
123 declareProperty("keepUnmatch", m_keepUnmatch= 0);
124 declareProperty("salvageTrk", m_salvageTrk = false);
125 declareProperty("dropMultiHotInLayer",m_dropMultiHotInLayer = false);
126 declareProperty("dropTrkPt", m_dropTrkPt = -999.);
127 declareProperty("d0Cut", m_d0Cut = 999.);
128 declareProperty("z0Cut", m_z0Cut = 999.);
129
130 declareProperty("minMdcDigi", m_minMdcDigi = 0);
131 declareProperty("countPropTime",m_countPropTime = true);
132 declareProperty("addHitCut", m_addHitCut = 5.);
133 declareProperty("dropHitsSigma",m_dropHitsSigma);
134 declareProperty("helixFitCut", m_helixFitCut);
135 declareProperty("minTrkProb", m_minTrkProb = 0.01);
136 declareProperty("csmax4", m_csmax4 = 50.);
137 declareProperty("csmax3", m_csmax3 = 1.);
138 declareProperty("helixFitSigma",m_helixFitSigma= 5.);
139 declareProperty("maxRcsInAddSeg",m_maxRcsInAddSeg= 50.);
140 declareProperty("nSigAddHitTrk", m_nSigAddHitTrk= 5.);
141 declareProperty("maxProca", m_maxProca= 0.6);
142 declareProperty("doSag", m_doSag= false);
143 declareProperty("lineFit", m_lineFit = false);
144 //declareProperty("cosmicFit", m_cosmicFit= false);
145}
146
147//--------------
148// Destructor --
149//--------------
151 delete m_bfield;
152}
153
155 // Get Mdc Detector Geometry
156 m_gm = MdcDetector::instance(m_doSag);
157 if(NULL == m_gm) return StatusCode::FAILURE;
159
160 return StatusCode::SUCCESS;
161}
162//--------------
163// Operations --
164//--------------
166 MsgStream log(msgSvc(), name());
167 log << MSG::INFO << "in initialize()" << endreq;
168
169 t_nTkTot = 0;
170 for (int i=0;i<20;i++) t_nTkNum[i]=0;
171
172 //m_flags.readPar(m_paramFile);
173#ifdef MDCXTIMEDEBUG
174 StatusCode tsc = service( "BesTimerSvc", m_timersvc);
175 if( tsc.isFailure() ) {
176 log << MSG::WARNING << name() << ": Unable to locate BesTimer Service" << endreq;
177 return StatusCode::FAILURE;
178 }
179 m_timer[0] = m_timersvc->addItem("Execution");
180 m_timer[0]->propName("Execution");
181 m_timer[1] = m_timersvc->addItem("findSeg");
182 m_timer[1]->propName("findSeg");
183 m_timer[2] = m_timersvc->addItem("findTrack");
184 m_timer[2]->propName("findTrack");
185 m_timer[3] = m_timersvc->addItem("fitting");
186 m_timer[3]->propName("fitting");
187#endif
188
189 if (m_helixFitCut.size() == 43){
190 for(int i=0; i<43; i++){
191 //MdcTrkReconCut_helix_fit[i] = m_helixFitCut[i];
192 TrkHelixFitter::nSigmaCut[i] = m_helixFitCut[i];
193 }
194 }
195 MdcxParameters::debug = m_debug;
196 MdcxParameters::minTrkProb = m_minTrkProb;
197 MdcxParameters::csmax4 = m_csmax4;
198 MdcxParameters::csmax3 = m_csmax3;
199 MdcxParameters::helixFitSigma = m_helixFitSigma;
200 MdcxParameters::maxRcsInAddSeg= m_maxRcsInAddSeg;
201 MdcxParameters::nSigAddHitTrk = m_nSigAddHitTrk;
202 MdcxParameters::maxProca = m_maxProca;
203 TrkHelixFitter::m_debug = (m_debug>7);
204 Pdt::readMCppTable(m_pdtFile);
205 MdcxFittedHel::debug = m_debug;
206
207
208 // Get MdcCalibFunSvc
209 IMdcCalibFunSvc* imdcCalibSvc;
210 StatusCode sc = service ("MdcCalibFunSvc", imdcCalibSvc);
211 m_mdcCalibFunSvc = dynamic_cast<MdcCalibFunSvc*>(imdcCalibSvc);
212 if ( sc.isFailure() ){
213 log << MSG::FATAL << "Could not load MdcCalibFunSvc!" << endreq;
214 return StatusCode::FAILURE;
215 }
216 MdcxHit::setMdcCalibFunSvc(m_mdcCalibFunSvc);
217 MdcxHit::setCountPropTime(m_countPropTime);
218
219 // Get RawDataProviderSvc
220 IRawDataProviderSvc* iRawDataProvider;
221 sc = service ("RawDataProviderSvc", iRawDataProvider);
222 if ( sc.isFailure() ){
223 log << MSG::FATAL << "Could not load RawDataProviderSvc!" << endreq;
224 return StatusCode::FAILURE;
225 }
226 m_rawDataProviderSvc = dynamic_cast<RawDataProviderSvc*>(iRawDataProvider);
227
228
229 //Initailize magnetic filed
230 sc = service ("MagneticFieldSvc",m_pIMF);
231 if(sc!=StatusCode::SUCCESS) {
232 log << MSG::ERROR << "Unable to open Magnetic field service"<<endreq;
233 }
234 m_bfield = new BField(m_pIMF);
235 log << MSG::INFO << "field z = "<<m_bfield->bFieldNominal()<< endreq;
236 m_context = new TrkContextEv(m_bfield);
237
238 if (m_hist) {bookNTuple();}
239 if (m_dropHitsSigma.size()==43){
240 for (int ii=0;ii<43;ii++) {
241 MdcxParameters::dropHitsSigma[ii]=m_dropHitsSigma[ii];
242 }
243 }
244
245 return StatusCode::SUCCESS;
246}
247
249
250
251 b_saveEvent=false;
252 setFilterPassed(b_saveEvent);
253#ifdef MDCXTIMEDEBUG
254 m_timer[0]->start();
255 m_timer[1]->start();
256#endif
257 MsgStream log(msgSvc(), name());
258 log << MSG::INFO << "in execute()" << endreq;
259 StatusCode sc;
260
261 nTk = 0; t_nTdsTk = 0; t_nDigi=0; t_nSeg=0; //yzhang for fill
262 //------------------------------------
263 // Get event No.
264 //------------------------------------
265 SmartDataPtr<Event::EventHeader> evtHead(eventSvc(),"/Event/EventHeader");
266 if (!evtHead) {
267 log << MSG::FATAL<< "Could not retrieve event header" << endreq;
268 return StatusCode::FAILURE;
269 }
270 m_eventNo = evtHead->eventNumber();
271 if(m_debug>0)std::cout << "x evt: "<<evtHead->runNumber()<<" " << m_eventNo<< std::endl;
272 long t_evtNo = m_eventNo;
273 g_eventNo = m_eventNo;
274 //if (t_evtNo % 1000 == 0) std::cout << "x evt: " << t_evtNo << std::endl;
275 IDataManagerSvc *dataManSvc;
276 DataObject *aTrackCol;
277 DataObject *aHitCol;
278 if (!m_salvageTrk) {
279 SmartIF<IDataManagerSvc> dataManSvc(eventSvc());
280 eventSvc()->findObject("/Event/Recon/RecMdcTrackCol",aTrackCol);
281 if(aTrackCol != NULL) {
282 dataManSvc->clearSubTree("/Event/Recon/RecMdcTrackCol");
283 eventSvc()->unregisterObject("/Event/Recon/RecMdcTrackCol");
284 }
285 eventSvc()->findObject("/Event/Recon/RecMdcHitCol",aHitCol);
286 if(aHitCol != NULL) {
287 dataManSvc->clearSubTree("/Event/Recon/RecMdcHitCol");
288 eventSvc()->unregisterObject("/Event/Recon/RecMdcHitCol");
289 }
290 }
291
292 //------------------------------------
293 // Initialize track collection in TDS
294 //------------------------------------
295 DataObject *aReconEvent;
296 eventSvc()->findObject("/Event/Recon",aReconEvent);
297 if (aReconEvent==NULL) {
298 aReconEvent = new ReconEvent();
299 sc = eventSvc()->registerObject("/Event/Recon",aReconEvent);
300 if(sc != StatusCode::SUCCESS) {
301 log << MSG::FATAL << "Could not register ReconEvent" <<endreq;
302 return StatusCode::FAILURE;
303 }
304 }
305 RecMdcTrackCol* trackList;
306 eventSvc()->findObject("/Event/Recon/RecMdcTrackCol",aTrackCol);
307 if (aTrackCol) {
308 trackList = dynamic_cast<RecMdcTrackCol*> (aTrackCol);
309 }else{
310 trackList = new RecMdcTrackCol;
311 sc = eventSvc()->registerObject(EventModel::Recon::RecMdcTrackCol, trackList);
312 if(!sc.isSuccess()) {
313 log << MSG::FATAL << " Could not register RecMdcTrack collection" <<endreq;
314 return StatusCode::FAILURE;
315 }
316 }
317 RecMdcHitCol* hitList;
318 eventSvc()->findObject("/Event/Recon/RecMdcHitCol",aHitCol);
319 if (aHitCol) {
320 hitList = dynamic_cast<RecMdcHitCol*> (aHitCol);
321 }else{
322 hitList = new RecMdcHitCol;
323 sc = eventSvc()->registerObject(EventModel::Recon::RecMdcHitCol, hitList);
324 if(!sc.isSuccess()) {
325 log << MSG::FATAL << " Could not register RecMdcHit collection" <<endreq;
326 return StatusCode::FAILURE;
327 }
328 }
329
330 //------------------------------------
331 // Initialize hit collection in TDS
332 //------------------------------------
333 DataObject *pnode = 0;
334 sc = eventSvc()->retrieveObject("/Event/Hit",pnode);
335 if(!sc.isSuccess()) {
336 pnode = new DataObject;
337 sc = eventSvc()->registerObject("/Event/Hit",pnode);
338 if(!sc.isSuccess()) {
339 log << MSG::FATAL << " Could not register /Event/Hit branch " <<endreq;
340 return StatusCode::FAILURE;
341 }
342 }
343 SmartDataPtr<MdcHitCol> m_hitCol(eventSvc(), "/Event/Hit/MdcHitCol");
344 if (!m_hitCol){
345 m_hitCol= new MdcHitCol;
346 sc = eventSvc()->registerObject("/Event/Hit/MdcHitCol",m_hitCol);
347 if(!sc.isSuccess()) {
348 log << MSG::FATAL << " Could not register hit collection" <<endreq;
349 return StatusCode::FAILURE;
350 }
351 }
352
353
354 //------------------------------------
355 // Get bunch time t0 (ns) and timing
356 //------------------------------------
357 m_bunchT0 = -999.;
358 SmartDataPtr<RecEsTimeCol> aevtimeCol(eventSvc(),"/Event/Recon/RecEsTimeCol");
359 if (!aevtimeCol || aevtimeCol->size()==0) {
360 log << MSG::WARNING<< "evt "<<m_eventNo<<" Could not find RecEsTimeCol"<< endreq;
361 return StatusCode::SUCCESS;
362 }
363
364 RecEsTimeCol::iterator iter_evt = aevtimeCol->begin();
365 for(; iter_evt!=aevtimeCol->end(); iter_evt++){
366 m_bunchT0 = (*iter_evt)->getTest();
367 m_t0Stat = (*iter_evt)->getStat();
368 if(m_debug>1) std::cout<<name()<<" "<<t_evtNo<<" t0 "<<m_bunchT0<<" t0Stat "<<m_t0Stat<<std::endl;
369 if ((m_t0Stat==0) || (m_bunchT0 < 0.) || (m_bunchT0 > 9999.0) ){
370 log << MSG::WARNING << "Skip evt:"<<m_eventNo<< " by t0 = "<<m_bunchT0 << endreq;
371 //return StatusCode::SUCCESS;
372 }
373 }
374 if(m_debug>1) std::cout<<name()<<" evt "<<t_evtNo<<" t0 "<<m_bunchT0<<" t0Stat "<<m_t0Stat<<std::endl;
375 int trigtiming=-10;
376 SmartDataPtr<TrigData> trigData(eventSvc(),"/Event/Trig/TrigData");
377 if(trigData){
378 log << MSG::INFO <<"Trigger conditions 0--43:"<<endreq;
379 for(int i = 0; i < 48; i++) {
380 log << MSG::INFO << trigData->getTrigCondName(i)<<" ---- "<<trigData->getTrigCondition(i)<< endreq;
381 }
382 for(int i = 0; i < 16; i++) log << MSG::INFO << "Trigger channel "<< i << ": " << trigData->getTrigChannel(i) << endreq;
383 m_timing=trigData->getTimingType();
384 //cout<<"-----------------trigger timing type-----------------------: "<<trigtiming<<endl;
385 log << MSG::INFO <<"Tigger Timing type: "<< trigtiming << endreq;
386 }
387
388 //------------------------------------
389 // Initialize MdcxHits
390 //------------------------------------
391 m_mdcxHits.reset();
392 uint32_t getDigiFlag = 0;
393 getDigiFlag += m_maxMdcDigi;
394 if(m_dropHot||m_salvageTrk) getDigiFlag |= MdcRawDataProvider::b_dropHot;
395 if(m_keepBadTdc) getDigiFlag |= MdcRawDataProvider::b_keepBadTdc;
396 if(m_keepUnmatch) getDigiFlag |= MdcRawDataProvider::b_keepUnmatch;
397 mdcDigiVec = m_rawDataProviderSvc->getMdcDigiVec(getDigiFlag);
398 t_nDigi = mdcDigiVec.size();
399
400 // fill Mc truth
401 //if(m_hist) fillMcTruth();
402
403 //skip event by hit numbe
404 if (t_nDigi < m_minMdcDigi){
405 if (0 == t_nDigi ){
406 log << MSG::WARNING << " No hits in MdcDigiVec" << endreq;
407 }
408 log << MSG::WARNING << " Skip this event for MdcDigiVec.size() < "<<m_minMdcDigi << endreq;
409 return StatusCode::SUCCESS;
410 }
411 m_mdcxHits.create(mdcDigiVec, m_bunchT0, m_cresol);
412 const HepAList<MdcxHit>& dchitlist = m_mdcxHits.GetMdcxHitList();
413
414 if(m_debug>2) m_mdcxHits.print(std::cout,6796);
415
416 //--------------------------------------------
417 // Make segments (MdcxSeg's) out of MdcxHit's
418 //--------------------------------------------
419 MdcxFindSegs dcsegs(dchitlist,m_debug);
420 const HepAList<MdcxSeg>& seglist = dcsegs.GetMdcxSeglist();
421 if(m_debug > 1 ){ dumpMdcxSegs(seglist);}
422 t_nSeg = seglist.length();
423 //if(m_hist){ fillMdcxSegs(seglist);}
424
425#ifdef MDCXTIMEDEBUG
426 m_timer[1]->stop();
427 m_timer[2]->start();
428#endif
429 //--------------------------------------------
430 // Make tracks (MdcxFittedHel's) out of MdcxSeg's
431 //--------------------------------------------
432 MdcxFindTracks dctrks(seglist,m_debug);
434 if(m_debug>1) dumpTrackList(firsttrkl);
435
436#ifdef MDCXTIMEDEBUG
437 m_timer[2]->stop();
438 m_timer[3]->start();
439#endif
440 //if(m_hist){ fillTrkl(firsttrkl);}
441
442 //if(m_debug>1){
443 // std::cout << "dchitlist after find tracks before MergeDups, nhits=" << dchitlist.length() << std::endl;
444 // for (int ii = 0; ii < dchitlist.length(); ii++) {
445 // dchitlist[ii]->print(std::cout, ii);
446 // }
447 // std::cout<<std::endl;
448 //}
449 MdcxMergeDups dcmergeem(firsttrkl,m_debug);
451
452 //if (m_debug > 1 ){
453 // cout << "MdcxTrackFinder: after MergeDups, have "
454 // << trkl.length() << " track(s). nhits=" << dchitlist.length() << endl;
455 // for (int ii = 0; ii < dchitlist.length(); ii++) {
456 // dchitlist[ii]->print(std::cout, ii);
457 // }
458 // std::cout<<std::endl;
459 //}
460
461 //---------------------------------------------------------
462 // Put my tracks into official fitter and store to TDS
463 //----------------------------------------------------------
464
465
466 sc = FitMdcxTrack(trkl, dchitlist, m_hitCol, trackList, hitList);
467 if (!sc.isSuccess()) {return StatusCode::SUCCESS;}
468 t_nTdsTk = trackList->size();
469
470 t_nTkTot += trackList->size();
471 if(t_nTdsTk<20) t_nTkNum[t_nTdsTk]++;
472
473#ifdef MDCXTIMEDEBUG
474 m_timer[0]->stop();
475 m_timer[3]->stop();
476#endif
477
478 if(m_hist) fillEvent();
479 if (m_debug > 0) {
480 DataObject* pNode;
481 eventSvc()->retrieveObject("/Event/Recon/RecMdcTrackCol",pNode);
482 RecMdcTrackCol *tmpTrackCol = dynamic_cast<RecMdcTrackCol*> (pNode);
483 eventSvc()->retrieveObject("/Event/Recon/RecMdcHitCol",pNode);
484 int nTdsTk = 0;
485 if(tmpTrackCol) nTdsTk = tmpTrackCol->size();
486
487 //if (t_evtNo % 1000 == 0) {
488 std::cout<< "MdcxTrackFinder: evtNo "<< m_eventNo << " t0="<<m_bunchT0
489 <<" Found " <<trkl.length()
490 <<" keep "<< t_nTdsTk
491 <<" finialy keep "<< nTdsTk;
492
493 int ndelete =0; trkl.length() - trackList->size();
494 if( ndelete>0 ) std::cout <<" delete "<< ndelete;
495 std::cout <<" track(s)" <<endl;
496 //}
497
498 if(m_debug>1)dumpTdsTrack(tmpTrackCol);
499 if(m_debug>1)dumpTrack(tmpTrackCol);
500 //dumpTdsHits(tmpHitCol);
501
502 }
503 if((trackList->size()!=4) ) b_saveEvent = true;
504 setFilterPassed(b_saveEvent);
505 return StatusCode::SUCCESS;
506}
507
509 MsgStream log(msgSvc(), name());
510 log << MSG::INFO << "in finalize()" << endreq;
511
512 std::cout<< " --after " << name() << " keep " << t_nTkTot<< " tracks "<<std::endl;
513 for(int i=0; i<20; i++){
514 if(t_nTkNum[i]>0)std::cout<< " nTk="<< i << " "<< t_nTkNum[i]<< std::endl;
515 }
516
517 //tot evtNo, trkNum
518 return StatusCode::SUCCESS;
519}
520
521
522void MdcxTrackFinder::MdcxHitsToHots(MdcxHel& mdcxHelix, const HepAList<MdcxHit>& mdcxHits,
523 TrkHitList* m_trkHitList, MdcHitCol* mdcHitCol) {
524
525 if ( 0 == mdcxHits.length() ) return;
526
527 int ihits = 0;
528 while(mdcxHits[ihits]) {
529 int ambig = 0;//=mdcxHelix.Doca_Samb();//yzhang delete 2011-10-13
530 const MdcHit* newhit = mdcxHits[ihits]->getMdcHit();
531 if ( 0 == newhit ) {
532 const MdcDigi* theDigi = mdcxHits[ihits]->getDigi();
533 int layer = MdcID::layer(mdcxHits[ihits]->getDigi()->identify());
534 int wire = MdcID::wire(mdcxHits[ihits]->getDigi()->identify());
535 m_digiMap[layer][wire] = mdcxHits[ihits]->getDigi();
536 MdcHit *thehit = new MdcHit(theDigi, m_gm);
537 thehit->setCalibSvc(m_mdcCalibFunSvc);
538 thehit->setCountPropTime(m_countPropTime);
539 //thehit->setCosmicFit(m_cosmicFit);
540
541 mdcHitCol->push_back(thehit);
542 newhit = thehit;
543 }
544 MdcRecoHitOnTrack temp(*newhit, ambig, m_bunchT0);//m_bunchT0 nano second here
545 MdcHitOnTrack* newhot = &temp;
546 //double fltLen = mdcxHelix.Doca_FLen();
547 //std::cout<<" fltLen-- "<<fltLen<<" "<<std::endl;
548 //newhot->setFltLen(fltLen);
549 //newhot->print(std::cout);
550 //std::cout<< __FILE__ << " ambig " << ambig << " append "<<std::endl;
551
552 // Store MdcxHits to TrkHitList
553 m_trkHitList->appendHot(newhot); //yzhang TEMP FIXME
554 newhot->setUsability(true);//yzhang add 2011-10-12 TEMP
555 newhot->setActivity(true);//yzhang add 2011-10-12 TEMP
556
557 ihits++;
558 }
559}
560
561StatusCode MdcxTrackFinder::FitMdcxTrack(HepAList<MdcxFittedHel>& trkl,
562 const HepAList<MdcxHit>& dchitlist, MdcHitCol* m_hitCol,
563 RecMdcTrackCol* trackList, RecMdcHitCol* hitList){
564 StatusCode sc;
565 MsgStream log(msgSvc(), name());
566
567 //--Add Hit to MdcxFittedHel
568 //if(m_debug>1){
569 // std::cout << "dchitlist before addHits nhits=" << dchitlist.length() << std::endl;
570 // for (int ii = 0; ii < dchitlist.length(); ii++) {
571 // dchitlist[ii]->print(std::cout, ii);
572 // }
573 // std::cout<<std::endl;
574 //}
575 MdcxAddHits dcaddem(trkl, dchitlist, m_addHitCut);
576 //if (m_debug > 1){
577 // cout << "MdcxTrackFinder: after AddHits, have "
578 // << trkl.length() << " track(s). nhits=" << dchitlist.length() << endl;
579 // for (int ii = 0; ii < dchitlist.length(); ii++) {
580 // dchitlist[ii]->print(std::cout, ii);
581 // }
582 // std::cout<<std::endl;
583 //}
584
585 TrkLineMaker linefactory;
586 TrkHelixMaker helixfactory;
587
588
589 int tkLen = trkl.length();//FIXME
590 for (int kk=0; kk< tkLen; kk++){
591 const HepAList<MdcxHit>& xhits = trkl[kk]->XHitList();
592 if(m_debug>2){
593 std::cout<<__FILE__<<" FitMdcxTrack "<<kk<< std::endl;
594 for(int i=0; i<xhits.length(); i++){ xhits[i]->print(std::cout); }
595 std::cout<<std::endl;
596 }
597
598 if(m_debug>2) std::cout<<" before add hits nhits="<<xhits.length()<< std::endl;
599 HepAList<MdcxHit> xass = dcaddem.GetAssociates(kk);
600 if(m_debug>2){
601 std::cout<<" after,add "<<xass.length()<<" hits"<<std::endl;
602 }
603
604 MdcxFittedHel mdcxHelix = *trkl[kk];
605 double thechisq = mdcxHelix.Chisq();
606 TrkExchangePar tt(-mdcxHelix.D0(),mdcxHelix.Phi0(),mdcxHelix.Omega(),-mdcxHelix.Z0(),-mdcxHelix.Tanl());
607 TrkRecoTrk* aTrack;
608 int nparm;
609
610 if (m_lineFit) {
611 /// m_bunchT0 in "ns" here, but "second" in TrkLineMaker&TrkHelixMaker
612 aTrack = linefactory.makeTrack(tt, thechisq, *m_context, m_bunchT0*1.e-9);
613 nparm = 4;
614 linefactory.setFlipAndDrop(*aTrack, true, true);
615 } else {
616 aTrack = helixfactory.makeTrack(tt, thechisq, *m_context, m_bunchT0*1.e-9);
617 nparm = 5;
618 helixfactory.setFlipAndDrop(*aTrack, true, true);
619 }
620
621 TrkHitList* m_trkHitList = aTrack->hits();
622 if (0 == m_trkHitList) {
623 delete aTrack;
624 aTrack = NULL;
625 continue;
626 }
627
628
629 MdcxHitsToHots(mdcxHelix, xhits, m_trkHitList, m_hitCol);
630 //std::cout<<"xhits---------------------"<<std::endl;//debug
631 //m_trkHitList->hotList().printAll(cout);//debug
632 MdcxHitsToHots(mdcxHelix, xass, m_trkHitList, m_hitCol);
633 //std::cout<<"xass----------------------"<<std::endl;//debug
634 //m_trkHitList->hotList().printAll(cout);//debug
635 //std::cout<<"size "<<m_trkHitList->hotList().nHit()<<std::endl;
636 //int beforDrop = m_trkHitList->hotList().nHit();
637 if(m_dropMultiHotInLayer) dropMultiHotInLayer(m_trkHitList);//yzhang debug FIXME
638 //int afterDrop = m_trkHitList->hotList().nHit();
639 //std::cout<<"drop "<<beforDrop-afterDrop<<" keep:"<<afterDrop<<::endl;
640 if(m_debug>1){ std::cout<< "== put to official fitter "<<std::endl;}
641 TrkErrCode err = m_trkHitList->fit();
642
643 if(m_debug>1){
644 std::cout<< "== after official fitter "<<std::endl;
645 aTrack->printAll(std::cout);
646 }
647
648 const TrkFit* theFit = aTrack->fitResult();
649 float rcs = 10000.0;
650
651 if (theFit) {
652 int ndof = theFit->nActive()-nparm;
653 if (ndof > 0) rcs = theFit->chisq()/float(ndof);
654 if (m_debug>1) {
655 if (4 == nparm) cout << " TrkLineMaker";
656 else cout << " TrkHelixMaker";
657 cout << " trkNo. " << kk << " success " << err.success() << " rcs " << rcs
658 << " chi2 " << theFit->chisq() << " nactive " << theFit->nActive() << endl;
659 }
660 }
661
662 if ( (1 == err.success()) && (rcs < 20.0) ) {
663 if(m_debug>1) std::cout<<"== fit success "<<std::endl;//yzhang debug
664 if (4 == nparm) {
665 linefactory.setFlipAndDrop(*aTrack, false, false);
666 } else {
667 helixfactory.setFlipAndDrop(*aTrack, false, false);
668 }
669 //-------------Stick the found tracks into the list in RecMdcTrackCol--------
670 if (m_debug>1) { cout << "MdcxTrackFinder: accept a track " << endl; }
671 // update history
672 aTrack->status()->addHistory(err,name().c_str()); //yzhang FIXME
673 store(aTrack, trackList, hitList);//aTrack have been deleted
674 } else if ( (2 == err.success()) && (rcs < 150.0) ) { //////// zoujh
675 if(m_debug > 1) std::cout<<"== fit success = 2, refit now"<<std::endl;//yzhang debug
676 int nrefit = 0;
677 while (nrefit++ < 5) {
678 if (m_debug>1) std::cout << "refit time " << nrefit << std::endl;
679 err = m_trkHitList->fit();
680 if (err.success() == 1) break;
681 }
682 if (err.success() == 1) {
683 if (4 == nparm) {
684 linefactory.setFlipAndDrop(*aTrack, false, false);
685 } else {
686 helixfactory.setFlipAndDrop(*aTrack, false, false);
687 }
688 //-------------Stick the found tracks into the list in RecMdcTrackCol--------
689 //if (m_debug>1) { cout << "MdcxTrackFinder: accept a track and store to TDS" << endl; }
690 // update history
691 aTrack->status()->addHistory(err,name().c_str()); //yzhang FIXME
692 store(aTrack, trackList, hitList);//aTrack have been deleted
693 } //////////////////////////zoujh
694 } else {
695 if (m_debug >1) {
696 std::cout<<"== fit faild "<<std::endl;//yzhang debug
697 err.print(cout);
698 cout << endl;
699 }
700 delete aTrack;
701 aTrack = NULL;
702 //---------------------------------------
703 // Fit no good; try a better input helix
704 //---------------------------------------
705 if(m_debug>1) std::cout<<"== Fit no good; try a better input helix"<<std::endl;//yzhang debug
706 mdcxHelix.Grow(*trkl[kk],xass);
707 mdcxHelix.VaryRes();
708 mdcxHelix.SetChiDofBail(1500);//yzhang add 2009-11-03
709 int fail = mdcxHelix.ReFit();
710 if(m_debug>1)std::cout<<__FILE__<<" refit fail:"<<fail<< std::endl;
711 if (!mdcxHelix.Fail()) {
712 const HepAList<MdcxHit>& bxhits = mdcxHelix.XHitList();
713 thechisq = mdcxHelix.Chisq();
714 TrkExchangePar tb(mdcxHelix.D0(),mdcxHelix.Phi0(),mdcxHelix.Omega(),mdcxHelix.Z0(),mdcxHelix.Tanl());
715 TrkRecoTrk* bTrack;
716 if (4 == nparm){
717 bTrack = linefactory.makeTrack(tb,thechisq,*m_context,m_bunchT0*1.e-9);
718 linefactory.setFlipAndDrop(*bTrack, false, false);
719 }else{
720 bTrack = helixfactory.makeTrack(tb,thechisq,*m_context,m_bunchT0*1.e-9);
721 helixfactory.setFlipAndDrop(*bTrack, false, false);
722 }
723 TrkHitList* bhits = bTrack->hits();
724 if (0 == bhits){delete bTrack; bTrack = NULL; continue;}
725
726 MdcxHitsToHots(mdcxHelix, bxhits, bhits, m_hitCol);
727 TrkErrCode berr = bhits->fit();
728 const TrkFit* bFit = bTrack->fitResult();
729 rcs=10000.0;
730 if (bFit) {
731 int ndof = bFit->nActive() - nparm;
732 if (ndof > 0) rcs = bFit->chisq()/float(ndof);
733 if (m_debug >1) {
734 if (4 == nparm) cout << " TrkLineMaker";
735 else cout << " TrkHelixMaker";
736 cout << " success trkNo. " << kk << " status " << berr.success() << " rcs " << rcs
737 << " chi2 " << bFit->chisq() << " nactive "<< bFit->nActive() << endl;
738 }
739 }
740 if ( ( 1 == berr.success() ) && ( rcs < 50.0 ) ) {
741 // update history
742 bTrack->status()->addHistory(berr,name().c_str());//yzhang FIXME
743 if (m_debug>1) {
744 cout << "MdcxTrackFinder: accept b track and store to TDS" << endl;
745 bTrack->printAll(cout);
746 }
747 store(bTrack, trackList, hitList);//bTrack have been deleted
748 } else {
749 if (m_debug>1) {
750 cout<< " fit failed "<<endl;
751 berr.print(cout);
752 cout << endl;
753 }
754 if (bTrack!=NULL) { delete bTrack; bTrack = NULL; }
755 }
756 }else{
757 //cout<< " grow and refit failed "<<endl;
758 }
759 }
760 }
761 if(m_debug >1) dumpTdsTrack(trackList);
762 return StatusCode::SUCCESS;
763
764}// end of FitMdcxTrack
765
766void MdcxTrackFinder::store(TrkRecoTrk* aTrack, RecMdcTrackCol *trackList,
767 RecMdcHitCol *hitList) {
768 assert (aTrack != NULL);
769 nTk++;
770 int trackId = trackList->size();
771 TrkExchangePar helix = aTrack->fitResult()->helix(0.);
772 if(m_dropTrkPt>0. && (aTrack->fitResult()->pt()<m_dropTrkPt)) {
773 if(m_debug>1) std::cout<<"MdcxTrackFinder delete track by pt "
774 <<aTrack->fitResult()->pt()<<"<ptCut "<<m_dropTrkPt << std::endl;
775 return;
776 }
777
778 if( ( (fabs(helix.d0())>m_d0Cut) ||( fabs(helix.z0())>m_z0Cut) ) ){
779 if(m_debug>1) std::cout<< name()<<" delete track by d0 "<<helix.d0()<<">d0Cut "<<m_d0Cut
780 <<" or z0 "<<helix.z0()<<" >z0Cut "<<m_z0Cut << std::endl;
781 return;
782 }
783
784 if(m_hist) fillTrack(aTrack);
785 MdcTrack mdcTrack(aTrack);//aTrack have been deleted in ~MdcTrack()
786 //tkStat: 0,PatRec 1,MdcxReco 2,Tsf 3,CurlFinder -1,Combined cosmic
787 int tkStat = 1;
788 int nHitbefore = hitList->size();
789
790 mdcTrack.storeTrack(trackId, trackList, hitList, tkStat);
791 int nHitAfter = hitList->size();
792 if (nHitAfter - nHitbefore <10 ) setFilterPassed(true);
793}
794
795void MdcxTrackFinder::dumpTrack(RecMdcTrackCol* trackList){
796 RecMdcTrackCol::iterator i_tk = trackList->begin();
797 for (;i_tk != trackList->end(); i_tk++) {
798 printTrack(*(i_tk));
799 }
800}// dumpTrack
801
802void MdcxTrackFinder::printTrack(RecMdcTrack* tk){
803 //yzhang debug
804 std::cout<< " MdcTrack Id:"<<tk->trackId() <<" q:"<< tk->charge()<< std::endl;
805 std::cout<< "dr Fi0 Cpa Dz Tanl Chi2 Ndf nSt FiTerm poca" << std::endl;
806 std::cout<<"(" <<setw(5) << tk->helix(0)<<","<< setw(5) << tk->helix(1)<<"," << setw(5) << tk->helix(2) <<","
807 << setw(5) << tk->helix(3) << ","<< setw(5) << tk->helix(4) <<")"
808 << setw(5) << tk->chi2() << setw(4) << tk->ndof()
809 << setw(4) << tk->getNhits() << setw(4) << tk->nster()
810 << setw(5) << tk->getFiTerm() <<tk->poca()<<std::endl;
811 std::cout<< " ErrMat "<<tk->err() << std::endl;
812
813 std::cout<< "hitId tkId (l,w) fltLen lr dt ddl tdc "//<<chi2Add
814 <<"doca entr z tprop stat " << std::endl;
815
816 HitRefVec hl = tk->getVecHits();
817 HitRefVec::iterator it = hl.begin();
818 for (;it!=hl.end();++it){
819 RecMdcHit* h = *it;
820 int layer = MdcID::layer(h->getMdcId());
821 double _vprop = (layer<8) ? Constants::vpropInner : Constants::vpropOuter;
822 const MdcLayer* _layerPtr = m_gm->Layer(layer);
823 double _zlen = _layerPtr->zLength();
824 double z = h->getZhit();
825 double tprop;
826 if (0 == layer%2){
827 tprop = (0.5*_zlen + z)/_vprop; //odd
828 }else{
829 tprop = (0.5*_zlen - z)/_vprop; //even
830 }
831 // build the sense wires
832 //const MdcSWire* wire = m_gm->Wire(MdcID::layer(h->getMdcId()),MdcID::wire(h->getMdcId()));
833 //double z = wire->zForw();
834 //while(z<wire->zRear()){
835 //HepPoint3D pos;
836 //Hep3Vector dir;
837 //if(!(wire->getTraj()==NULL)){
838 //wire->getTraj()->getInfo(z,pos,dir);
839 //}
840 //std::cout<<"("<< wire->layer()->layNum()<<","<<wire->cell()<<" "<<wire->Id()<<")";
841 //std::cout<<" z, sag:"<<z
842 //<<", "<<wire->getTraj()->deltaY(z-wire->zForw())<<std::endl;
843 ////<<" pos:"<<pos <<" dir:"<<dir<<std::endl;
844 //z+=1.;
845 //}
846 std::cout<< setw(3) << h->getId() << setw(2) << h->getTrkId() <<
847 "(" << MdcID::layer(h->getMdcId()) <<"," << MdcID::wire(h->getMdcId()) <<")"<<
848 setw(10) << h->getFltLen() <<
849 setw(2) << h->getFlagLR() <<setw(10) << h->getDriftT() <<
850 setw(12) << h->getDriftDistLeft() <<//setw(8) << h->getErrDriftDistRight() <<
851 setw(8) << h->getTdc() <<//setw(10) << h->getChisqAdd() <<
852 setw(10) << h->getDoca() <<setw(10) << h->getEntra() <<
853 setw(10) << h->getZhit() << setw(10) << tprop<<
854 setw(2)<< h->getStat() << std::endl;
855 }
856}//print track
857
858void MdcxTrackFinder::bookNTuple(){
859 MsgStream log(msgSvc(), name());
860 StatusCode sc;
861 if(!m_hist) return;
862
863 g_poison = histoSvc()->book( "poison", "poison",43,0,42,288,0,288 );
864 g_csmax4 = histoSvc()->book( "csmax4", "csmax4",100,0,100 );
865 g_csmax3 = histoSvc()->book( "csmax3", "csmax3",50000,0,500000 );
866 g_omegag = histoSvc()->book( "omegag", "omegag",1000 ,0.,1.);
867 g_dPhiAU = histoSvc()->book( "dPhiAU", "dPhiAU",1000 ,0.,3.5);
868 g_dPhiAU_0 = histoSvc()->book( "dPhiAU_0", "dPhiAU_0",1000 ,0.,3.5);
869 g_dPhiAU_1 = histoSvc()->book( "dPhiAU_1", "dPhiAU_1",1000 ,0.,3.5);
870 g_dPhiAU_5 = histoSvc()->book( "dPhiAU_5", "dPhiAU_5",1000 ,0.,3.5);
871 g_dPhiAU_7 = histoSvc()->book( "dPhiAU_7", "dPhiAU_7",1000 ,0.,3.5);
872 g_dPhiAV = histoSvc()->book( "dPhiAV", "dPhiAV",1000 ,0.,3.5);
873 g_addSegPhi = histoSvc()->book( "addSegPhi", "addSegPhi",1000 ,0.,3.5);
874 g_dPhiAV_0 = histoSvc()->book( "dPhiAV_0", "dPhiAV_0",1000 ,0.,3.5);
875 g_dPhiAV_1 = histoSvc()->book( "dPhiAV_1", "dPhiAV_1",1000 ,0.,3.5);
876 g_dPhiAV_6 = histoSvc()->book( "dPhiAV_6", "dPhiAV_6",1000 ,0.,3.5);
877 g_dPhiAV_8 = histoSvc()->book( "dPhiAV_8", "dPhiAV_8",1000 ,0.,3.5);
878 g_trkllmk = histoSvc()->book( "trkllmk", "trkllmk",43,0.,42);
879 g_trklgood = histoSvc()->book( "trklgood", "trklgood",43,0.,42);
880 g_trklcircle = histoSvc()->book( "trklcircle", "trklcircle",43,0.,42);
881 g_trklhelix= histoSvc()->book( "trklhelix", "trklhelix",43,0.,42);
882 g_trkldrop1= histoSvc()->book( "trkldrop1", "trkldrop1",43,0.,42);
883 g_trkldrop2= histoSvc()->book( "trkldrop2", "trkldrop2",43,0.,42);
884 g_trklappend1= histoSvc()->book( "trklappend1", "trklappend1",43,0.,42);
885 g_trklappend2= histoSvc()->book( "trklappend2", "trklappend2",43,0.,42);
886 g_trklappend3= histoSvc()->book( "trklappend3", "trklappend3",43,0.,42);
887 //g_fitOmega= histoSvc()->book( "fitOmega", "fitOmega",200,0.,100.);
888 g_trklfirstProb= histoSvc()->book( "trklfirstProb", "trklfirstProb",200,0.,2.);
889 g_trkltemp= histoSvc()->book( "trkltemp", "trkltemp",200,-3.5,3.5);
890
891 g_trklproca= histoSvc()->book( "trklproca", "trklproca",100,0.,5.);
892 g_trkle= histoSvc()->book( "trkle", "trkle",200,0.,0.025);
893 g_trkld= histoSvc()->book( "trkld", "trkld",200,-1.2,1.2);
894 g_trklel= histoSvc()->book( "trklel", "trklel",200,0.,0.025,43,0,43);
895 g_trkldl= histoSvc()->book( "trkldl", "trkldl",200,-1.2,1.2,43,0,43);
896 g_trkldoca= histoSvc()->book( "trkldoca", "trkldoca",200,-1.2,1.2);
897 g_trkllayer= histoSvc()->book( "trkllayer", "trkllayer",43,0,43);
898 g_dropHitsSigma= histoSvc()->book( "dropHitsSigma", "dropHitsSigma",43,0,43,100,0,11);
899 g_addHitCut= histoSvc()->book( "addHitCut", "addHitCut",100,0,10);
900 g_addHitCut2d= histoSvc()->book( "addHitCut2d", "addHitCut2d",43,0,43,100,0,10);
901 //g_addSegPhiDiff = histoSvc()->book( "addSegPhiDiff", "addSegPhiDiff",100,-6.28,6.28);
902
903 NTuplePtr nt1(ntupleSvc(), "MdcxReco/rec");
904 if ( nt1 ) { m_xtuple1 = nt1;}
905 else {
906 m_xtuple1 = ntupleSvc()->book ("MdcxReco/rec", CLID_ColumnWiseTuple, "MdcxReco reconsturction results");
907 if ( m_xtuple1 ) {
908 sc = m_xtuple1->addItem ("t0", m_xt0);
909 sc = m_xtuple1->addItem ("timing", m_xtiming);
910 sc = m_xtuple1->addItem ("t0Stat", m_xt0Stat);
911 sc = m_xtuple1->addItem ("t0Truth", m_xt0Truth);
912
913 sc = m_xtuple1->addItem ("nSlay", m_xnSlay);
914 sc = m_xtuple1->addItem ("p", m_xp);
915 sc = m_xtuple1->addItem ("pt", m_xpt);
916 sc = m_xtuple1->addItem ("pz", m_xpz);
917 sc = m_xtuple1->addItem ("d0", m_xd0);
918 sc = m_xtuple1->addItem ("phi0", m_xphi0);
919 sc = m_xtuple1->addItem ("cpa", m_xcpa);
920 sc = m_xtuple1->addItem ("z0", m_xz0);
921 sc = m_xtuple1->addItem ("tanl", m_xtanl);
922 sc = m_xtuple1->addItem ("q", m_xq);
923 sc = m_xtuple1->addItem ("pocax", m_xpocax);
924 sc = m_xtuple1->addItem ("pocay", m_xpocay);
925 sc = m_xtuple1->addItem ("pocaz", m_xpocaz);
926
927 sc = m_xtuple1->addItem ("evtNo", m_xevtNo);
928 sc = m_xtuple1->addItem ("nSt", m_xnSt);
929 sc = m_xtuple1->addItem ("nAct", m_xnAct);
930 sc = m_xtuple1->addItem ("nDof", m_xnDof);
931 sc = m_xtuple1->addItem ("chi2", m_xchi2);
932 sc = m_xtuple1->addItem ("tkId", m_xtkId);
933 sc = m_xtuple1->addItem ("nHit", m_xnHit, 0, 7000);//# of hit/rec track
934 sc = m_xtuple1->addIndexedItem ("resid", m_xnHit, m_xresid);
935 sc = m_xtuple1->addIndexedItem ("sigma", m_xnHit, m_xsigma);
936 sc = m_xtuple1->addIndexedItem ("driftD", m_xnHit, m_xdriftD);
937 sc = m_xtuple1->addIndexedItem ("driftT", m_xnHit, m_xdriftT);
938 sc = m_xtuple1->addIndexedItem ("doca", m_xnHit, m_xdoca);
939 sc = m_xtuple1->addIndexedItem ("entra", m_xnHit, m_xentra);
940 sc = m_xtuple1->addIndexedItem ("ambig", m_xnHit, m_xambig);
941 sc = m_xtuple1->addIndexedItem ("fltLen", m_xnHit, m_xfltLen);
942 sc = m_xtuple1->addIndexedItem ("tof", m_xnHit, m_xtof);
943 sc = m_xtuple1->addIndexedItem ("act", m_xnHit, m_xact);
944 sc = m_xtuple1->addIndexedItem ("tdc", m_xnHit, m_xtdc);
945 sc = m_xtuple1->addIndexedItem ("adc", m_xnHit, m_xadc);
946 sc = m_xtuple1->addIndexedItem ("layer", m_xnHit, m_xlayer);
947 sc = m_xtuple1->addIndexedItem ("wire", m_xnHit, m_xwire);
948 sc = m_xtuple1->addIndexedItem ("x", m_xnHit, m_xx);
949 sc = m_xtuple1->addIndexedItem ("y", m_xnHit, m_xy);
950 sc = m_xtuple1->addIndexedItem ("z", m_xnHit, m_xz);
951 } else { // did not manage to book the N tuple....
952 log << MSG::ERROR << " Cannot book N-tuple: MdcxReco/rec" << endmsg;
953 //return;
954 }
955 }//end book of nt1
956 NTuplePtr nt4(ntupleSvc(), "MdcxReco/evt");
957 if ( nt4 ) { m_xtupleEvt = nt4;}
958 else {
959 m_xtupleEvt = ntupleSvc()->book ("MdcxReco/evt", CLID_ColumnWiseTuple, "MdcxReco event data");
960 if ( m_xtupleEvt ) {
961 sc = m_xtupleEvt->addItem ("evtNo", m_xt4EvtNo );
962 sc = m_xtupleEvt->addItem ("nRecTk", m_xt4nRecTk );
963 sc = m_xtupleEvt->addItem ("nTdsTk", m_xt4nTdsTk );
964 sc = m_xtupleEvt->addItem ("t0", m_xt4t0);
965 sc = m_xtupleEvt->addItem ("t0Stat", m_xt4t0Stat);
966 sc = m_xtupleEvt->addItem ("t0Truth", m_xt4t0Truth);
967 sc = m_xtupleEvt->addItem ("time", m_xt4time);
968 sc = m_xtupleEvt->addItem ("timeSeg", m_xt4timeSeg);
969 sc = m_xtupleEvt->addItem ("timeTrack", m_xt4timeTrack);
970 sc = m_xtupleEvt->addItem ("timeFit", m_xt4timeFit);
971 sc = m_xtupleEvt->addItem ("nSeg", m_xt4nSeg);
972 sc = m_xtupleEvt->addItem ("nDigi", m_xt4nDigi, 0, 7000);//# of hit/mc track
973 sc = m_xtupleEvt->addIndexedItem ("layer", m_xt4nDigi, m_xt4Layer);
974 sc = m_xtupleEvt->addIndexedItem ("rt", m_xt4nDigi, m_xt4Time);
975 sc = m_xtupleEvt->addIndexedItem ("rc", m_xt4nDigi, m_xt4Charge);
976 //sc = m_xtupleEvt->addIndexedItem ("rawHit", m_xt4nDigi, m_xt4rawHit);
977 //sc = m_xtupleEvt->addIndexedItem ("recHit", m_xt4nDigi, m_xt4recHit);
978 } else { // did not manage to book the N tuple....
979 log << MSG::ERROR << "Cannot book N-tuple: MdcxReco/evt" << endmsg;
980 //return;
981 }
982 }// end of book nt4
983
984 //book tuple of segment
985 NTuplePtr ntSeg(ntupleSvc(), "MdcxReco/seg");
986 if ( ntSeg ) { m_xtupleSeg = ntSeg;}
987 else {
988 m_xtupleSeg = ntupleSvc()->book ("MdcxReco/seg", CLID_ColumnWiseTuple, "MdcxTrackFinder segment data");
989 if ( m_xtupleSeg ) {
990 sc = m_xtupleSeg->addItem ("sl", m_xtsSl);
991 sc = m_xtupleSeg->addItem ("d0", m_xtsD0);
992 sc = m_xtupleSeg->addItem ("omega", m_xtsOmega);
993 sc = m_xtupleSeg->addItem ("phi0", m_xtsPhi0);
994 sc = m_xtupleSeg->addItem ("d0sl", m_xtsD0_sl_approx);
995 sc = m_xtupleSeg->addItem ("phi0sl", m_xtsPhi0_sl_approx);
996 sc = m_xtupleSeg->addItem ("xbbrrf", m_xtsXline_bbrrf);
997 sc = m_xtupleSeg->addItem ("ybbrrf", m_xtsYline_bbrrf);
998 sc = m_xtupleSeg->addItem ("xslope", m_xtsXline_slope);
999 sc = m_xtupleSeg->addItem ("yslope", m_xtsYline_slope);
1000 sc = m_xtupleSeg->addItem ("pat", m_xtsPat);
1001 sc = m_xtupleSeg->addItem ("chisq", m_xtsChisq);
1002 sc = m_xtupleSeg->addItem ("nDigi", m_xtsNDigi, 0, 20);
1003 sc = m_xtupleSeg->addIndexedItem ("layer", m_xtsNDigi, m_xtsLayer);
1004 sc = m_xtupleSeg->addIndexedItem ("wire", m_xtsNDigi, m_xtsWire);
1005 sc = m_xtupleSeg->addIndexedItem ("inSeg", m_xtsNDigi, m_xtsInSeg);
1006 }else{
1007 log << MSG::ERROR << "Cannot book N-tuple: MdcxReco/evt" << endmsg;
1008 //return;
1009 }
1010 }
1011 NTuplePtr nt5(ntupleSvc(), "MdcxReco/trkl");
1012 if ( nt5 ) { m_xtupleTrkl= nt5;}
1013 else {
1014 m_xtupleTrkl= ntupleSvc()->book ("MdcxReco/trkl", CLID_RowWiseTuple, "MdcxReco track info");
1015 if ( m_xtupleTrkl) {
1016 sc = m_xtupleTrkl->addItem ("layer", m_xt5Layer);
1017 sc = m_xtupleTrkl->addItem ("wire", m_xt5Wire);
1018 }else{
1019 log << MSG::ERROR << "Cannot book N-tuple: MdcxReco/trkl" << endmsg;
1020 //return;
1021 }
1022 }
1023
1024 NTuplePtr ntCsmcSew(ntupleSvc(), "MdcxReco/csmc");
1025 if ( ntCsmcSew ) { m_xtupleCsmcSew = ntCsmcSew;}
1026 else {
1027 m_xtupleCsmcSew = ntupleSvc()->book ("MdcxReco/csmc", CLID_ColumnWiseTuple, "MdcxReco reconsturction results");
1028 if ( m_xtupleCsmcSew ) {
1029 sc = m_xtupleCsmcSew->addItem ("dD0", m_csmcD0);
1030 sc = m_xtupleCsmcSew->addItem ("dPhi0", m_csmcPhi0);
1031 sc = m_xtupleCsmcSew->addItem ("dZ0", m_csmcZ0);
1032 sc = m_xtupleCsmcSew->addItem ("dOmega", m_csmcOmega);
1033 sc = m_xtupleCsmcSew->addItem ("dPt", m_csmcPt);
1034 sc = m_xtupleCsmcSew->addItem ("dTanl", m_csmcTanl);
1035 }else{
1036 log << MSG::ERROR << "Cannot book N-tuple: MdcxReco/csmc" << endmsg;
1037 //return;
1038 }
1039 }
1040 NTuplePtr ntSegAdd1(ntupleSvc(), "MdcxReco/addSeg1");
1041 if ( ntSegAdd1 ) { m_xtupleAddSeg1 = ntSegAdd1;}
1042 else {
1043 m_xtupleAddSeg1 = ntupleSvc()->book ("MdcxReco/addSeg1", CLID_ColumnWiseTuple, "MdcxReco event data");
1044 if ( m_xtupleAddSeg1 ) {
1045 sc = m_xtupleAddSeg1->addItem ("same", m_addSegSame);
1046 sc = m_xtupleAddSeg1->addItem ("seedSl", m_addSegSeedSl );
1047 sc = m_xtupleAddSeg1->addItem ("seedPhi", m_addSegSeedPhi );
1048 sc = m_xtupleAddSeg1->addItem ("seedPhiLay",m_addSegSeedPhiLay );
1049 sc = m_xtupleAddSeg1->addItem ("seedLen", m_addSegSeedLen );
1050 sc = m_xtupleAddSeg1->addItem ("seedD0", m_addSegSeedD0 );
1051 sc = m_xtupleAddSeg1->addItem ("seedPhi0", m_addSegSeedPhi0 );
1052 sc = m_xtupleAddSeg1->addItem ("addSl", m_addSegAddSl );
1053 sc = m_xtupleAddSeg1->addItem ("addPhi", m_addSegAddPhi );
1054 sc = m_xtupleAddSeg1->addItem ("addPhiLay", m_addSegAddPhiLay );
1055 sc = m_xtupleAddSeg1->addItem ("addLen", m_addSegAddLen );
1056 sc = m_xtupleAddSeg1->addItem ("addD0", m_addSegAddD0 );
1057 sc = m_xtupleAddSeg1->addItem ("addPhi0", m_addSegAddPhi0 );
1058 }
1059 else{
1060 log << MSG::ERROR << "Cannot book N-tuple: MdcxReco/addSeg1" << endmsg;
1061 //return;
1062 }
1063 }
1064
1065 NTuplePtr ntSegComb(ntupleSvc(), "MdcxReco/segComb");
1066 if ( ntSegComb ) { m_xtupleSegComb = ntSegComb;}
1067 else {
1068 m_xtupleSegComb = ntupleSvc()->book ("MdcxReco/segComb", CLID_ColumnWiseTuple, "MdcxReco event data");
1069 if ( m_xtupleSegComb ) {
1070 sc = m_xtupleSegComb->addItem ("evtNo", m_segCombEvtNo);
1071 sc = m_xtupleSegComb->addItem ("omega", m_segCombOmega);
1072 sc = m_xtupleSegComb->addItem ("sameAU", m_segCombSameAU);
1073 sc = m_xtupleSegComb->addItem ("sameUV", m_segCombSameUV);
1074 sc = m_xtupleSegComb->addItem ("dlenAU", m_segCombDLenAU);
1075 sc = m_xtupleSegComb->addItem ("dlenUV", m_segCombDLenUV);
1076 sc = m_xtupleSegComb->addItem ("slA", m_segCombSlA);
1077 sc = m_xtupleSegComb->addItem ("slU", m_segCombSlU);
1078 sc = m_xtupleSegComb->addItem ("slV", m_segCombSlV);
1079 sc = m_xtupleSegComb->addItem ("phiA", m_segCombPhiA);
1080 sc = m_xtupleSegComb->addItem ("phiU", m_segCombPhiU);
1081 sc = m_xtupleSegComb->addItem ("phiV", m_segCombPhiV);
1082 }
1083 else{
1084 log << MSG::ERROR << "Cannot book N-tuple: MdcxReco/segComb" << endmsg;
1085 //return;
1086 }
1087 }
1088
1089 NTuplePtr ntDropHits(ntupleSvc(), "MdcxReco/dropHits");
1090 if ( ntDropHits ) { m_xtupleDropHits = ntDropHits;}
1091 else {
1092 m_xtupleDropHits = ntupleSvc()->book ("MdcxReco/dropHits", CLID_ColumnWiseTuple, "MdcxReco event data");
1093 if ( m_xtupleDropHits ) {
1094 sc = m_xtupleDropHits->addItem ("evtNo", m_segDropHitsEvtNo);
1095 sc = m_xtupleDropHits->addItem ("layer", m_segDropHitsLayer);
1096 sc = m_xtupleDropHits->addItem ("wire", m_segDropHitsWire);
1097 sc = m_xtupleDropHits->addItem ("pull", m_segDropHitsPull);
1098 sc = m_xtupleDropHits->addItem ("doca", m_segDropHitsDoca);
1099 sc = m_xtupleDropHits->addItem ("sigma", m_segDropHitsSigma);
1100 sc = m_xtupleDropHits->addItem ("drift", m_segDropHitsDrift);
1101 sc = m_xtupleDropHits->addItem ("mcTkId", m_segDropHitsMcTkId);
1102 } else{
1103 log << MSG::ERROR << "Cannot book N-tuple: MdcxReco/dropHits" << endmsg;
1104 //return;
1105 }
1106 }
1107 NTuplePtr ntAddSeg2(ntupleSvc(), "MdcxReco/addSeg2");
1108 if ( ntAddSeg2 ) { m_xtupleAddSeg2 = ntAddSeg2;}
1109 else {
1110 m_xtupleAddSeg2 = ntupleSvc()->book ("MdcxReco/addSeg2", CLID_ColumnWiseTuple, "MdcxReco event data");
1111 if ( m_xtupleAddSeg2 ) {
1112 sc = m_xtupleAddSeg2->addItem ("evtNo", m_addSegEvtNo);
1113 sc = m_xtupleAddSeg2->addItem ("slayer", m_addSegSlayer);
1114 sc = m_xtupleAddSeg2->addItem ("poca", m_addSegPoca);
1115 sc = m_xtupleAddSeg2->addItem ("len", m_addSegLen);
1116 } else{
1117 log << MSG::ERROR << "Cannot book N-tuple: MdcxReco/addSeg2" << endmsg;
1118 //return;
1119 }
1120 }
1121 //--------------end of book ntuple------------------
1122}
1123
1124
1125void MdcxTrackFinder::fillEvent(){
1126 //-----------get raw digi-----------------------
1127 if (m_xtupleEvt == NULL ) return;
1128#ifdef MDCXTIMEDEBUG
1129 m_xt4time = m_timer[0]->elapsed();
1130 m_xt4timeSeg = m_timer[1]->elapsed();
1131 m_xt4timeTrack = m_timer[2]->elapsed();
1132 m_xt4timeFit = m_timer[3]->elapsed();
1133#endif
1134 m_xt4EvtNo = m_eventNo;
1135 m_xt4t0 = m_bunchT0;
1136 m_xt4t0Stat = m_t0Stat;
1137 m_xt4t0Truth = m_t0Truth;
1138 m_xt4nRecTk = nTk;
1139 m_xt4nTdsTk = t_nTdsTk;
1140 m_xt4nDigi = t_nDigi;
1141 m_xt4nSeg = t_nSeg;
1142 int iDigi=0;
1143 MdcDigiCol::iterator iter = mdcDigiVec.begin();
1144 for (;iDigi<t_nDigi; iter++ ) {
1145 int l = MdcID::layer((*iter)->identify());
1146 m_xt4Layer[iDigi] = l;
1147 //int w = MdcID::wire((*iter)->identify());
1148 m_xt4Time[iDigi] = RawDataUtil::MdcTime((*iter)->getTimeChannel());
1149 m_xt4Charge[iDigi] = RawDataUtil::MdcCharge((*iter)->getChargeChannel());
1150 //m_xt4rawHit[l]++;
1151 iDigi++;
1152 }//end for iter
1153 m_xtupleEvt->write();
1154}
1155
1156void MdcxTrackFinder::fillTrack(TrkRecoTrk* atrack){
1157
1158 if(!m_xtuple1) return;
1159 //-----------get MC truth data-------------------
1160 m_xevtNo = m_eventNo;
1161 int recHitMap[43][288]={0};
1162 //int nTk = (*m_tracks).nTrack();
1163 //for (int itrack = 0; itrack < nTk; itrack++) {
1164 if (atrack==NULL) return;
1165
1166 const TrkFit* fitResult = atrack->fitResult();
1167 if (fitResult == 0) return;//check the fit worked
1168
1169 //for fill ntuples
1170 int nSt=0; //int nSeg=0;
1171 int seg[11] = {0}; int segme;
1172 //-----------hit list-------------
1173 HitRefVec hitRefVec;
1174 const TrkHitList* hitlist = atrack->hits();
1175
1176 TrkHitList::hot_iterator hot = hitlist->begin();
1177 int layerCount[43]={0};
1178 int i=0;
1179 for (;hot!= hitlist->end();hot++){
1180
1181 const MdcRecoHitOnTrack* recoHot;
1182 recoHot = dynamic_cast<const MdcRecoHitOnTrack*> (&(*hot));
1183 int layer = recoHot->mdcHit()->layernumber();
1184 int wire = recoHot->mdcHit()->wirenumber();
1185 m_xlayer[i] = layer;
1186 //m_xt4recHit[layer]++;//fill rec hit for hit eff
1187 m_xwire[i] = wire;
1188 layerCount[layer]++;
1189 recHitMap[layer][wire]++;
1190 m_xambig[i] = recoHot->wireAmbig();// wire ambig
1191 //fltLen
1192 double fltLen = recoHot->fltLen();
1193 m_xfltLen[i] = fltLen;
1194 double tof = recoHot->getParentRep()->arrivalTime(fltLen);
1195 //position
1196 HepPoint3D pos; Hep3Vector dir;
1197 recoHot->getParentRep()->traj().getInfo(fltLen,pos,dir);
1198 m_xx[i] = pos.x();
1199 m_xy[i] = pos.y();
1200 m_xz[i] = pos.z();
1201 m_xdriftT[i] = recoHot->mdcHit()->driftTime(tof,pos.z());
1202 m_xtof[i] = tof;
1203 m_xdriftD[i] = recoHot->drift();
1204 m_xsigma[i] = recoHot->hitRms();
1205 m_xtdc[i] = recoHot->rawTime();
1206 m_xadc[i] = recoHot->mdcHit()->charge();
1207 m_xdoca[i] = recoHot->dcaToWire();//sign w.r.t. dirft() FIXME
1208 m_xentra[i] = recoHot->entranceAngle();
1209 //m_xentraHit[i] = recoHot->entranceAngleHit();
1210 m_xact[i] = recoHot->isActive();
1211 //resid
1212 double res=999.,rese=999.;
1213 if (recoHot->resid(res,rese,false)){
1214 }else{}
1215 m_xresid[i] = res;
1216 //for n seg
1217 segme=0;
1218 if ( layer >0 ) {segme=(layer-1)/4;}
1219 seg[segme]++;
1220 if (recoHot->layer()->view()) { ++nSt; }
1221 i++;
1222 }// end fill hit
1223
1224 int nSlay=0;
1225 for(int i=0;i<11;i++){
1226 if (seg[i]>0) nSlay++;
1227 }
1228 m_xnSlay = nSlay;
1229
1230 //------------fill track result-------------
1231 //m_xtkId = itrack;
1232 //track parameters at closest approach to beamline
1233 double fltLenPoca = 0.0;
1234 TrkExchangePar helix = fitResult->helix(fltLenPoca);
1235 m_xphi0 = helix.phi0();
1236 m_xtanl = helix.tanDip();
1237 m_xz0 = helix.z0();
1238 m_xd0 = helix.d0();
1239 if(fabs(m_xz0)>20||fabs(m_xd0)>2) {
1240 //b_saveEvent = true;
1241 if(m_debug>1) std::cout<<"evt:"<<m_eventNo<<" BigVtx:"
1242 <<" d0 "<<helix.d0()<<" z0 "<<helix.z0()<<std::endl;
1243 }
1244 m_xpt = fitResult->pt();
1245 if (m_xpt > 0.00001){
1246 m_xcpa = fitResult->charge()/fitResult->pt();
1247 }
1248 //momenta and position
1249 CLHEP::Hep3Vector p1 = fitResult->momentum(fltLenPoca);
1250 double px,py,pz,pxy;
1251 pxy = fitResult->pt();
1252 px = p1.x();
1253 py = p1.y();
1254 pz = p1.z();
1255 m_xp = p1.mag();
1256 m_xpz = pz;
1257 HepPoint3D poca = fitResult->position(fltLenPoca);
1258 m_xpocax = poca.x();
1259 m_xpocay = poca.y();
1260 m_xpocaz = poca.z();
1261
1262 m_xq = fitResult->charge();
1263 m_xnAct = fitResult->nActive();
1264 m_xnHit = hitlist->nHit();
1265 m_xnDof = fitResult->nDof();
1266 m_xnSt = nSt;
1267 m_xchi2 = fitResult->chisq();
1268 //for (int l=0;l<43;l++) m_xlayerCount[l] = layerCount[l];
1269 m_xt0 = m_bunchT0;
1270 m_xtiming = m_timing;
1271 m_xt0Stat = m_t0Stat;
1272 m_xt0Truth= m_t0Truth;
1273 m_xtuple1->write();
1274 //}//end of loop rec trk list
1275
1276}//end of
1277
1278void MdcxTrackFinder::dumpMdcxSegs(const HepAList<MdcxSeg>& segList) const {
1279 cout << name()<<" found " <<segList.length() << " segs :" << endl;
1280 for (int i =0; i< segList.length(); i++){
1281 std::cout<<i<<" ";
1282 segList[i]->printSeg();
1283 }
1284}
1285
1286void MdcxTrackFinder::fillMdcxSegs(const HepAList<MdcxSeg>& segList) const {
1287 if(!m_xtupleSeg) return;
1288
1289 int cellMax[43] ={
1290 40,44,48,56, 64,72,80,80, 76,76,88,88,
1291 100,100,112,112, 128,128,140,140, 160,160,160,160,
1292 176,176,176,176, 208,208,208,208, 240,240,240,240,
1293 256,256,256,256, 288,288,288 };
1294 // Fill hits of every layer after segment finding
1295 int hitInSegList[43][288];
1296 for (int ii=0;ii<43;ii++){
1297 for (int jj=0;jj<cellMax[ii];jj++){ hitInSegList[ii][jj] = 0; }
1298 }
1299 for (int i = 0; i < segList.length(); i++) {
1300 MdcxSeg* aSeg = segList[i];
1301 for (int ihit=0 ; ihit< aSeg->Nhits() ; ihit++){
1302 const MdcxHit* hit = aSeg->XHitList()[ihit];
1303 int layer = hit->Layer();
1304 int wire = hit->WireNo();
1305 hitInSegList[layer][wire]++;
1306 }
1307 }
1308 for (int i = 0; i < segList.length(); i++) {
1309 MdcxSeg* aSeg = segList[i];
1310 if(aSeg==NULL)continue;
1311 m_xtsSl = aSeg->SuperLayer();
1312 m_xtsD0 = aSeg->D0();
1313 m_xtsOmega = aSeg->Omega();
1314 m_xtsPhi0 = aSeg->Phi0();
1317 m_xtsXline_bbrrf = aSeg->Xline_bbrrf();
1318 m_xtsYline_bbrrf = aSeg->Yline_bbrrf();
1319 m_xtsXline_slope = aSeg->Xline_slope();
1320 m_xtsYline_slope = aSeg->Yline_slope();
1321 int patIndex = -1;
1322 for (int ii = 0;ii<14;ii++){
1323 if (aSeg->Pat() == (int) MdcxSegPatterns::patt4[ii]){
1324 patIndex = ii;
1325 break;
1326 }
1327 }
1328 for (int ii = 0;ii<20;ii++){
1329 if (aSeg->Pat() == (int) MdcxSegPatterns::patt3[ii]){
1330 patIndex = ii +15;
1331 break;
1332 }
1333 }
1334 m_xtsPat = patIndex;
1335 m_xtsChisq = aSeg->Chisq();
1336 m_xtsNDigi = aSeg->Nhits();
1337 for (int ihit=0 ; ihit< aSeg->Nhits() ; ihit++){
1338 const MdcxHit* hit = aSeg->XHitList()[ihit];
1339 int layer = hit->Layer();
1340 int wire = hit->WireNo();
1341 m_xtsLayer[ihit] = layer;
1342 m_xtsWire[ihit] = wire;
1343 m_xtsInSeg[ihit] = hitInSegList[layer][wire];
1344 //std::cout<< "hitInSegList "<<hitInSegList[layer][wire] << std::endl;//yzhang debug
1345 }
1346 m_xtupleSeg->write();
1347 }
1348
1349}//end of fillMdcxSegs
1350
1351void MdcxTrackFinder::dumpTrackList(const HepAList<MdcxFittedHel>& trackList) const {
1352 std::cout<< "dump track list " <<std::endl;
1353 for (int i=0;i<trackList.length();i++){
1354 std::cout<< "track "<<i<<std::endl;
1355 for (int ihit=0 ; ihit< trackList[i]->Nhits() ; ihit++){
1356 const MdcxHit* hit = trackList[i]->XHitList()[ihit];
1357 int layer = hit->Layer();
1358 int wire = hit->WireNo();
1359 std::cout<< " ("<<layer << ","<< wire<<") ";
1360 }
1361 std::cout<<std::endl;
1362 }
1363}// end of dumpTrackList
1364
1365void MdcxTrackFinder::fillTrkl(const HepAList<MdcxFittedHel>& firsttrkl) const {
1366 if(!m_xtupleTrkl) return;
1367 int nDigi = 0;
1368 int iDigi = 0;
1369 for (int i =0; i< firsttrkl.length(); i++){
1370 nDigi+= firsttrkl[i]->Nhits();
1371 }
1372 for (int i=0;i<firsttrkl.length();i++){
1373 for (int ihit=0 ; ihit< firsttrkl[i]->Nhits() ; ihit++){
1374 const MdcxHit* hit = firsttrkl[i]->XHitList()[ihit];
1375 int layer = hit->Layer();
1376 int wire = hit->WireNo();
1377 m_xt5Layer = layer;
1378 m_xt5Wire = wire;
1379 m_xtupleTrkl->write();
1380 iDigi++;
1381 }
1382 }
1383}//end of fillTrkl
1384
1385void MdcxTrackFinder::fillMcTruth(){
1386 MsgStream log(msgSvc(), name());
1387 StatusCode sc;
1388 //Initialize
1389 for (int ii=0;ii<43;ii++){
1390 for (int jj=0;jj<288;jj++){
1391 haveDigi[ii][jj] = -2;
1392 }
1393 }
1394 int nDigi = mdcDigiVec.size();
1395 for (int iDigi =0 ;iDigi<nDigi; iDigi++ ) {
1396 int l = MdcID::layer((mdcDigiVec[iDigi])->identify());
1397 int w = MdcID::wire((mdcDigiVec[iDigi])->identify());
1398 //haveDigi[l][w]=(mdcDigiVec[iDigi])->getTrackIndex();
1399 haveDigi[l][w]=1;
1400 }
1401
1402 if( m_mcHist ){
1403 //------------------get event start time truth-----------
1404 m_t0Truth = -10.;
1405 SmartDataPtr<Event::McParticleCol> mcParticleCol(eventSvc(),"/Event/MC/McParticleCol");
1406 if(!mcParticleCol){
1407 log << MSG::INFO << "Could not retrieve McParticelCol" << endreq;
1408 }else {
1409 Event::McParticleCol::iterator iter_mc = mcParticleCol->begin();
1410 for (;iter_mc != mcParticleCol->end(); iter_mc++){
1411 if ((*iter_mc)->primaryParticle()){
1412 m_t0Truth = (*iter_mc)->initialPosition().t();
1413 //px = (*iter_mc)->initialFourMomentum().x()/1000.;//GeV
1414 //py = (*iter_mc)->initialFourMomentum().y()/1000.;//GeV
1415 //pz = (*iter_mc)->initialFourMomentum().z()/1000.;//GeV
1416 }
1417 }
1418 }
1419 }
1420 //------------------Retrieve MC truth MdcMcHit------------
1421 /*
1422 SmartDataPtr<Event::MdcMcHitCol> mcMdcMcHitCol(eventSvc(),"/Event/MC/MdcMcHitCol");
1423 if (!mcMdcMcHitCol) {
1424 log << MSG::WARNING << "Could not find MdcMcHit" << endreq;
1425 }else{
1426 Event::MdcMcHitCol::iterator iter_mchit = mcMdcMcHitCol->begin();
1427 for (;iter_mchit != mcMdcMcHitCol->end(); iter_mchit++ ) {
1428 const Identifier id= (*iter_mchit)->identify();
1429 int layer = MdcID::layer(id);
1430 int wire = MdcID::wire(id);
1431 mcDrift[layer][wire] = (*iter_mchit)->getDriftDistance(); //drift in MC.
1432 mcLR[layer][wire] = (*iter_mchit)->getPositionFlag();
1433 mcX[layer][wire] = (*iter_mchit)->getPositionX();
1434 mcY[layer][wire] = (*iter_mchit)->getPositionY();
1435 mcZ[layer][wire] = (*iter_mchit)->getPositionZ();
1436 if (mcLR[layer][wire] == 0) mcLR[layer][wire] = -1;
1437 }
1438 }
1439 */
1440}
1441
1442void MdcxTrackFinder::dropMultiHotInLayer(TrkHitList* trkHitList){
1443 double tdr[43];
1444 double tdr_wire[43];
1445 for(int i=0; i<43; i++){tdr[i]=9999.;}
1446
1447 // make flag
1448 TrkHotList::hot_iterator hotIter = trkHitList->hotList().begin();
1449 while (hotIter!=trkHitList->hotList().end()) {
1450 MdcHitOnTrack* hot = const_cast<MdcHitOnTrack*> (&(*hotIter->mdcHitOnTrack()));
1451
1452 //driftTime(tof,z)
1453 double dt = hot->mdcHit()->driftTime(0.,0.);
1454 int layer = hot->mdcHit()->layernumber();
1455 int wire = hot->mdcHit()->wirenumber();
1456
1457 //std::cout<<__FILE__<<" "<<dt<< std::endl;
1458 if (dt < tdr[layer]) {
1459 tdr[layer] = dt;
1460 tdr_wire[layer] = wire;
1461 }
1462 hotIter++;
1463 }
1464
1465 //std::cout<<" tdr wire ";
1466 //for(int i=0;i<43;i++){
1467 //std::cout<<i<<","<<tdr_wire[i]<<" tdr="<<tdr[i]<<"\n";
1468 //}
1469 // inactive multi hit
1470 hotIter = trkHitList->hotList().begin();
1471 while (hotIter!=trkHitList->hotList().end()) {
1472 int layer = hotIter->mdcHitOnTrack()->mdcHit()->layernumber();
1473 int wire = hotIter->mdcHitOnTrack()->mdcHit()->wirenumber();
1474 //double dt = hotIter->mdcHitOnTrack()->mdcHit()->driftTime(0.,0.)-m_bunchT0;
1475
1476 if ((tdr[layer] <9998.) && (tdr_wire[layer]!=wire)){
1477 MdcHitOnTrack* hot = const_cast<MdcHitOnTrack*> (&(*hotIter->mdcHitOnTrack()));
1478 hot->setActivity(false);
1479 trkHitList->removeHit( hotIter->mdcHitOnTrack()->mdcHit() );
1480 //std::cout<<__FILE__<<" inactive "<< layer<<" "<<wire<<" dt "<<dt << std::endl;
1481 }else{
1482 hotIter++;
1483 }
1484 }
1485}
1486void MdcxTrackFinder::dumpTdsTrack(RecMdcTrackCol* trackList){
1487 std::cout<< "tksize = "<<trackList->size() << std::endl;//yzhang debug
1488 RecMdcTrackCol::iterator it = trackList->begin();
1489 for (;it!= trackList->end();it++){
1490 RecMdcTrack *tk = *it;
1491 std::cout<< "//====RecMdcTrack "<<tk->trackId()<<"====:" << std::endl;
1492 cout <<" d0 "<<tk->helix(0)
1493 <<" phi0 "<<tk->helix(1)
1494 <<" cpa "<<tk->helix(2)
1495 <<" z0 "<<tk->helix(3)
1496 <<" tanl "<<tk->helix(4)
1497 <<endl;
1498 std::cout<<" q "<<tk->charge()
1499 <<" theta "<<tk->theta()
1500 <<" phi "<<tk->phi()
1501 <<" x0 "<<tk->x()
1502 <<" y0 "<<tk->y()
1503 <<" z0 "<<tk->z()
1504 <<" r0 "<<tk->r()
1505 <<endl;
1506 std::cout <<" p "<<tk->p()
1507 <<" pt "<<tk->pxy()
1508 <<" px "<<tk->px()
1509 <<" py "<<tk->py()
1510 <<" pz "<<tk->pz()
1511 <<endl;
1512 std::cout<<" tkStat "<<tk->stat()
1513 <<" chi2 "<<tk->chi2()
1514 <<" ndof "<<tk->ndof()
1515 <<" nhit "<<tk->getNhits()
1516 <<" nst "<<tk->nster()
1517 <<endl;
1518
1519 int nhits = tk->getVecHits().size();
1520 std::cout<<nhits <<" Hits: " << std::endl;
1521 for(int ii=0; ii <nhits ; ii++){
1522 Identifier id(tk->getVecHits()[ii]->getMdcId());
1523 int layer = MdcID::layer(id);
1524 int wire = MdcID::wire(id);
1525 cout<<"("<< layer <<","<<wire<<","<<tk->getVecHits()[ii]->getStat()
1526 <<",lr:"<<tk->getVecHits()[ii]->getFlagLR()<<") ";
1527 }//end of hit list
1528 std::cout << " "<< std::endl;
1529 }//end of tk list
1530 std::cout << " "<< std::endl;
1531
1532
1533}//end of dumpTdsTrack
1534
1535void MdcxTrackFinder::dumpTdsHits(RecMdcHitCol* hitList){
1536 std::cout<<__FILE__<<" "<<__LINE__<<" All hits in TDS, nhit="<<hitList->size()<< std::endl;
1537 RecMdcHitCol::iterator it = hitList->begin();
1538 for(;it!= hitList->end(); it++){
1539 RecMdcHit* h = (*it);
1540 Identifier id(h->getMdcId());
1541 int layer = MdcID::layer(id);
1542 int wire = MdcID::wire(id);
1543 cout<<"("<< layer <<","<<wire<<") lr:"<<h->getFlagLR()<<" stat:"<<h->getStat()<<" tk:"<<h->getTrkId() <<" doca:"<<setw(10)<<h->getDoca()<<std::endl;
1544 }//end of hit list
1545}
double w
EvtStreamInputIterator< typename Generator::result_type > iter(Generator gen, int N=0)
TGraph2DErrors * dt
Definition: McCor.cxx:45
int recHitMap[43][288]
Definition: MdcHistItem.h:260
ObjectVector< MdcHit > MdcHitCol
Definition: MdcHit.h:129
NTuple::Array< double > m_xfltLen
Definition: MdcxHistItem.h:103
AIDA::IHistogram1D * g_dPhiAV
Definition: MdcxHistItem.h:30
NTuple::Item< double > m_addSegAddPhiLay
Definition: MdcxHistItem.h:170
NTuple::Array< double > m_xwire
Definition: MdcxHistItem.h:110
NTuple::Item< double > m_xtsXline_bbrrf
Definition: MdcxHistItem.h:139
NTuple::Item< double > m_xt0Stat
Definition: MdcxHistItem.h:70
NTuple::Array< double > m_xdriftT
Definition: MdcxHistItem.h:99
AIDA::IHistogram1D * g_trkldrop1
Definition: MdcxHistItem.h:40
NTuple::Array< double > m_xdriftD
Definition: MdcxHistItem.h:98
NTuple::Item< double > m_xpocax
Definition: MdcxHistItem.h:82
NTuple::Item< double > m_xtsChisq
Definition: MdcxHistItem.h:143
NTuple::Item< long > m_xnHit
Definition: MdcxHistItem.h:90
AIDA::IHistogram1D * g_dPhiAU_1
Definition: MdcxHistItem.h:27
NTuple::Item< double > m_xtsYline_slope
Definition: MdcxHistItem.h:142
NTuple::Item< long > m_xt5Layer
Definition: MdcxHistItem.h:149
NTuple::Item< double > m_xtsYline_bbrrf
Definition: MdcxHistItem.h:140
NTuple::Item< double > m_xq
Definition: MdcxHistItem.h:81
AIDA::IHistogram1D * g_csmax3
Definition: MdcxHistItem.h:23
AIDA::IHistogram1D * g_trklappend2
Definition: MdcxHistItem.h:43
NTuple::Item< double > m_addSegSeedPhi0
Definition: MdcxHistItem.h:167
NTuple::Array< double > m_xambig
Definition: MdcxHistItem.h:105
NTuple::Item< double > m_xtiming
Definition: MdcxHistItem.h:69
NTuple::Item< double > m_segCombDLenUV
Definition: MdcxHistItem.h:187
AIDA::IHistogram1D * g_trklappend1
Definition: MdcxHistItem.h:42
AIDA::IHistogram1D * g_dPhiAV_0
Definition: MdcxHistItem.h:31
NTuple::Item< long > m_addSegSlayer
Definition: MdcxHistItem.h:178
AIDA::IHistogram2D * g_poison
Definition: MdcxHistItem.h:21
NTuple::Item< double > m_xchi2
Definition: MdcxHistItem.h:88
NTuple::Array< double > m_xresid
Definition: MdcxHistItem.h:96
NTuple::Item< long > m_addSegSame
Definition: MdcxHistItem.h:161
NTuple::Item< double > m_xpocaz
Definition: MdcxHistItem.h:84
AIDA::IHistogram1D * g_trklproca
Definition: MdcxHistItem.h:49
NTuple::Item< long > m_segDropHitsWire
Definition: MdcxHistItem.h:198
NTuple::Item< long > m_xtsNDigi
Definition: MdcxHistItem.h:145
NTuple::Item< double > m_addSegSeedSl
Definition: MdcxHistItem.h:162
NTuple::Tuple * m_xtupleTrkl
Definition: MdcxHistItem.h:66
NTuple::Item< double > m_xnAct
Definition: MdcxHistItem.h:91
NTuple::Array< double > m_xtof
Definition: MdcxHistItem.h:104
NTuple::Tuple * m_xtupleAddSeg1
Definition: MdcxHistItem.h:160
NTuple::Item< double > m_addSegAddLen
Definition: MdcxHistItem.h:171
AIDA::IHistogram1D * g_addSegPhi
Definition: MdcxHistItem.h:35
AIDA::IHistogram1D * g_dPhiAV_1
Definition: MdcxHistItem.h:32
NTuple::Item< long > m_xtsSl
Definition: MdcxHistItem.h:133
NTuple::Item< double > m_xt4nTdsTk
Definition: MdcxHistItem.h:117
NTuple::Item< double > m_addSegLen
Definition: MdcxHistItem.h:179
AIDA::IHistogram1D * g_trkle
Definition: MdcxHistItem.h:52
NTuple::Item< long > m_xtsPat
Definition: MdcxHistItem.h:144
AIDA::IHistogram1D * g_dPhiAU_5
Definition: MdcxHistItem.h:28
NTuple::Item< double > m_xp
Definition: MdcxHistItem.h:72
AIDA::IHistogram1D * g_trkld
Definition: MdcxHistItem.h:51
NTuple::Array< long > m_xtsWire
Definition: MdcxHistItem.h:147
NTuple::Item< double > m_xt0
Definition: MdcxHistItem.h:68
AIDA::IHistogram1D * g_trkltemp
Definition: MdcxHistItem.h:47
NTuple::Item< double > m_segDropHitsDoca
Definition: MdcxHistItem.h:200
NTuple::Item< long > m_xt4nSeg
Definition: MdcxHistItem.h:122
NTuple::Item< double > m_xtsPhi0_sl_approx
Definition: MdcxHistItem.h:138
AIDA::IHistogram1D * g_csmax4
Definition: MdcxHistItem.h:22
NTuple::Item< double > m_csmcTanl
Definition: MdcxHistItem.h:158
NTuple::Item< double > m_segDropHitsSigma
Definition: MdcxHistItem.h:201
AIDA::IHistogram2D * g_trklel
Definition: MdcxHistItem.h:56
NTuple::Item< double > m_xtsD0
Definition: MdcxHistItem.h:134
NTuple::Item< double > m_addSegSeedPhiLay
Definition: MdcxHistItem.h:164
NTuple::Array< double > m_xact
Definition: MdcxHistItem.h:106
NTuple::Item< double > m_csmcOmega
Definition: MdcxHistItem.h:156
NTuple::Item< double > m_xpz
Definition: MdcxHistItem.h:74
NTuple::Item< double > m_addSegSeedLen
Definition: MdcxHistItem.h:165
NTuple::Item< double > m_xnSt
Definition: MdcxHistItem.h:86
NTuple::Array< long > m_xtsInSeg
Definition: MdcxHistItem.h:148
NTuple::Item< double > m_segCombSlV
Definition: MdcxHistItem.h:190
NTuple::Item< long > m_segDropHitsLayer
Definition: MdcxHistItem.h:197
NTuple::Item< double > m_addSegSeedD0
Definition: MdcxHistItem.h:166
AIDA::IHistogram1D * g_dPhiAU_7
Definition: MdcxHistItem.h:29
NTuple::Tuple * m_xtuple1
Definition: MdcxHistItem.h:63
AIDA::IHistogram2D * g_dropHitsSigma
Definition: MdcxHistItem.h:57
NTuple::Item< double > m_xtsPhi0
Definition: MdcxHistItem.h:136
AIDA::IHistogram1D * g_trklappend3
Definition: MdcxHistItem.h:44
NTuple::Item< double > m_xt4t0
Definition: MdcxHistItem.h:118
AIDA::IHistogram1D * g_trklhelix
Definition: MdcxHistItem.h:39
NTuple::Item< double > m_xt4timeFit
Definition: MdcxHistItem.h:126
NTuple::Item< double > m_xt4t0Truth
Definition: MdcxHistItem.h:120
NTuple::Item< double > m_segDropHitsMcTkId
Definition: MdcxHistItem.h:203
NTuple::Array< double > m_xdoca
Definition: MdcxHistItem.h:100
AIDA::IHistogram1D * g_trkllmk
Definition: MdcxHistItem.h:36
NTuple::Item< double > m_xevtNo
Definition: MdcxHistItem.h:85
AIDA::IHistogram1D * g_addHitCut
Definition: MdcxHistItem.h:58
AIDA::IHistogram2D * g_trkldl
Definition: MdcxHistItem.h:55
NTuple::Item< long > m_xt4t0Stat
Definition: MdcxHistItem.h:119
NTuple::Item< double > m_xt4nRecTk
Definition: MdcxHistItem.h:116
NTuple::Item< double > m_segCombPhiA
Definition: MdcxHistItem.h:191
NTuple::Item< long > m_xt4EvtNo
Definition: MdcxHistItem.h:115
AIDA::IHistogram1D * g_dPhiAU_0
Definition: MdcxHistItem.h:26
NTuple::Item< double > m_segCombSlA
Definition: MdcxHistItem.h:188
NTuple::Item< double > m_xt4timeTrack
Definition: MdcxHistItem.h:125
NTuple::Array< double > m_xsigma
Definition: MdcxHistItem.h:97
NTuple::Item< double > m_xt4time
Definition: MdcxHistItem.h:123
NTuple::Item< double > m_segCombDLenAU
Definition: MdcxHistItem.h:186
AIDA::IHistogram2D * g_addHitCut2d
Definition: MdcxHistItem.h:59
NTuple::Item< long > m_xt5Wire
Definition: MdcxHistItem.h:150
AIDA::IHistogram1D * g_trklgood
Definition: MdcxHistItem.h:38
AIDA::IHistogram1D * g_trkllayer
Definition: MdcxHistItem.h:54
NTuple::Item< double > m_xphi0
Definition: MdcxHistItem.h:77
NTuple::Item< double > m_segCombSlU
Definition: MdcxHistItem.h:189
NTuple::Item< long > m_segCombEvtNo
Definition: MdcxHistItem.h:182
NTuple::Item< double > m_segCombPhiU
Definition: MdcxHistItem.h:192
NTuple::Item< double > m_csmcPhi0
Definition: MdcxHistItem.h:154
NTuple::Array< double > m_xadc
Definition: MdcxHistItem.h:108
NTuple::Array< double > m_xt4Time
Definition: MdcxHistItem.h:128
AIDA::IHistogram1D * g_dPhiAU
Definition: MdcxHistItem.h:25
NTuple::Item< double > m_segCombSameUV
Definition: MdcxHistItem.h:185
NTuple::Item< double > m_addSegAddPhi
Definition: MdcxHistItem.h:169
NTuple::Item< double > m_xt4timeSeg
Definition: MdcxHistItem.h:124
NTuple::Item< double > m_addSegSeedPhi
Definition: MdcxHistItem.h:163
NTuple::Array< double > m_xx
Definition: MdcxHistItem.h:111
NTuple::Item< double > m_xz0
Definition: MdcxHistItem.h:79
NTuple::Item< double > m_csmcPt
Definition: MdcxHistItem.h:157
NTuple::Item< long > m_xnSlay
Definition: MdcxHistItem.h:76
AIDA::IHistogram1D * g_trkldrop2
Definition: MdcxHistItem.h:41
NTuple::Array< double > m_xtdc
Definition: MdcxHistItem.h:107
NTuple::Item< double > m_segCombPhiV
Definition: MdcxHistItem.h:193
NTuple::Item< double > m_xtanl
Definition: MdcxHistItem.h:80
NTuple::Item< double > m_addSegAddPhi0
Definition: MdcxHistItem.h:173
NTuple::Item< double > m_xt0Truth
Definition: MdcxHistItem.h:71
AIDA::IHistogram1D * g_dPhiAV_8
Definition: MdcxHistItem.h:34
NTuple::Tuple * m_xtupleSegComb
Definition: MdcxHistItem.h:181
NTuple::Item< double > m_addSegAddSl
Definition: MdcxHistItem.h:168
NTuple::Item< double > m_segCombOmega
Definition: MdcxHistItem.h:183
NTuple::Array< long > m_xt4Layer
Definition: MdcxHistItem.h:127
NTuple::Item< double > m_xtsD0_sl_approx
Definition: MdcxHistItem.h:137
NTuple::Item< double > m_xtkId
Definition: MdcxHistItem.h:89
NTuple::Item< double > m_addSegAddD0
Definition: MdcxHistItem.h:172
NTuple::Item< long > m_segDropHitsEvtNo
Definition: MdcxHistItem.h:196
NTuple::Item< double > m_segCombSameAU
Definition: MdcxHistItem.h:184
NTuple::Tuple * m_xtupleEvt
Definition: MdcxHistItem.h:64
NTuple::Item< double > m_xpt
Definition: MdcxHistItem.h:73
NTuple::Item< double > m_xpocay
Definition: MdcxHistItem.h:83
NTuple::Item< double > m_segDropHitsPull
Definition: MdcxHistItem.h:199
NTuple::Item< double > m_xcpa
Definition: MdcxHistItem.h:78
NTuple::Item< double > m_csmcD0
Definition: MdcxHistItem.h:153
AIDA::IHistogram1D * g_trklfirstProb
Definition: MdcxHistItem.h:46
NTuple::Tuple * m_xtupleAddSeg2
Definition: MdcxHistItem.h:175
AIDA::IHistogram1D * g_omegag
Definition: MdcxHistItem.h:24
NTuple::Array< long > m_xtsLayer
Definition: MdcxHistItem.h:146
AIDA::IHistogram1D * g_trkldoca
Definition: MdcxHistItem.h:53
NTuple::Array< double > m_xentra
Definition: MdcxHistItem.h:101
NTuple::Item< double > m_xnDof
Definition: MdcxHistItem.h:87
NTuple::Array< double > m_xlayer
Definition: MdcxHistItem.h:109
NTuple::Array< double > m_xy
Definition: MdcxHistItem.h:112
NTuple::Array< double > m_xz
Definition: MdcxHistItem.h:113
NTuple::Item< double > m_xtsOmega
Definition: MdcxHistItem.h:135
NTuple::Item< double > m_xd0
Definition: MdcxHistItem.h:75
NTuple::Item< long > m_addSegEvtNo
Definition: MdcxHistItem.h:176
NTuple::Item< double > m_segDropHitsDrift
Definition: MdcxHistItem.h:202
NTuple::Item< long > m_xt4nDigi
Definition: MdcxHistItem.h:121
AIDA::IHistogram1D * g_dPhiAV_6
Definition: MdcxHistItem.h:33
NTuple::Tuple * m_xtupleDropHits
Definition: MdcxHistItem.h:195
NTuple::Item< double > m_csmcZ0
Definition: MdcxHistItem.h:155
NTuple::Tuple * m_xtupleSeg
Definition: MdcxHistItem.h:65
NTuple::Item< double > m_addSegPoca
Definition: MdcxHistItem.h:177
NTuple::Tuple * m_xtupleCsmcSew
Definition: MdcxHistItem.h:152
AIDA::IHistogram1D * g_trklcircle
Definition: MdcxHistItem.h:37
NTuple::Item< double > m_xtsXline_slope
Definition: MdcxHistItem.h:141
NTuple::Array< double > m_xt4Charge
Definition: MdcxHistItem.h:129
double MdcTrkReconCut_helix_fit[43]
Definition: MdcTrkReconCut.h:6
ObjectVector< RecMdcHit > RecMdcHitCol
Definition: RecMdcHit.h:99
ObjectVector< RecMdcTrack > RecMdcTrackCol
Definition: RecMdcTrack.h:79
SmartRefVector< RecMdcHit > HitRefVec
Definition: RecMdcTrack.h:22
INTupleSvc * ntupleSvc()
IHistogramSvc * histoSvc()
IMessageSvc * msgSvc()
static const double vpropInner
Definition: Constants.h:49
static const double vpropOuter
Definition: Constants.h:50
const double theta() const
Definition: DstMdcTrack.h:59
const double r() const
Definition: DstMdcTrack.h:64
const double py() const
Definition: DstMdcTrack.h:56
const HepSymMatrix err() const
const double chi2() const
Definition: DstMdcTrack.h:66
const int charge() const
Definition: DstMdcTrack.h:53
const int trackId() const
Definition: DstMdcTrack.h:52
const double px() const
Definition: DstMdcTrack.h:55
const double phi() const
Definition: DstMdcTrack.h:60
const double pz() const
Definition: DstMdcTrack.h:57
const double pxy() const
Definition: DstMdcTrack.h:54
const int ndof() const
Definition: DstMdcTrack.h:67
const int stat() const
Definition: DstMdcTrack.h:65
const HepVector helix() const
......
const double z() const
Definition: DstMdcTrack.h:63
const double p() const
Definition: DstMdcTrack.h:58
const int nster() const
Definition: DstMdcTrack.h:68
const double y() const
Definition: DstMdcTrack.h:62
const HepPoint3D poca() const
Definition: DstMdcTrack.h:40
const double x() const
Definition: DstMdcTrack.h:61
const MdcLayer * Layer(unsigned id) const
Definition: MdcDetector.h:33
static MdcDetector * instance()
Definition: MdcDetector.cxx:21
virtual const MdcHit * mdcHit() const
int wireAmbig() const
double rawTime() const
double dcaToWire() const
double drift() const
Definition: MdcHitOnTrack.h:75
double entranceAngle() const
const MdcLayer * layer() const
Definition: MdcHit.h:44
void setCalibSvc(const MdcCalibFunSvc *calibSvc)
Definition: MdcHit.cxx:136
unsigned layernumber() const
Definition: MdcHit.h:61
unsigned wirenumber() const
Definition: MdcHit.h:62
double charge() const
Definition: MdcHit.h:65
double driftTime(double tof, double z) const
Definition: MdcHit.cxx:142
void setCountPropTime(const bool count)
Definition: MdcHit.h:86
static int layer(const Identifier &id)
Values of different levels (failure returns 0)
Definition: MdcID.cxx:49
static int wire(const Identifier &id)
Definition: MdcID.cxx:54
double zLength(void) const
Definition: MdcLayer.h:44
int view(void) const
Definition: MdcLayer.h:28
const MdcHit * mdcHit() const
const HepAList< MdcxSeg > & GetMdcxSeglist()
Definition: MdcxFindSegs.h:42
const HepAList< MdcxFittedHel > & GetMdcxTrklist()
int Nhits() const
Definition: MdcxFittedHel.h:49
const HepAList< MdcxHit > & XHitList() const
Definition: MdcxFittedHel.h:59
int SuperLayer(int hitno=0) const
static int debug
float Chisq() const
Definition: MdcxFittedHel.h:45
int Fail(float Probmin=0.0) const
void SetChiDofBail(float r)
MdcxFittedHel & Grow(const MdcxFittedHel &, const HepAList< MdcxHit > &)
double Phi0() const
Definition: MdcxHel.h:54
double D0() const
Definition: MdcxHel.h:53
double Omega() const
Definition: MdcxHel.h:55
double Z0() const
Definition: MdcxHel.h:56
double Tanl() const
Definition: MdcxHel.h:57
int Layer() const
Definition: MdcxHit.h:67
int WireNo() const
Definition: MdcxHit.h:66
static void setMdcDetector(const MdcDetector *gm)
Definition: MdcxHit.cxx:68
static void setMdcCalibFunSvc(const MdcCalibFunSvc *calibSvc)
Definition: MdcxHit.cxx:60
static void setCountPropTime(bool countPropTime)
Definition: MdcxHit.cxx:64
void print(std::ostream &o, int pmax=10) const
Definition: MdcxHits.cxx:72
const HepAList< MdcxHit > & GetMdcxHitList()
Definition: MdcxHits.h:40
void reset()
Definition: MdcxHits.cxx:55
void create(MdcDigiVec digiVec, float c0=0.0, float cresol=0.0180)
Definition: MdcxHits.cxx:59
const HepAList< MdcxFittedHel > & GetMergedTrklist()
Definition: MdcxMergeDups.h:33
static double maxProca
static int debug
static double maxRcsInAddSeg
static double minTrkProb
static float dropHitsSigma[43]
static double csmax4
static double helixFitSigma
static double csmax3
static double nSigAddHitTrk
static const unsigned patt4[14]
static const unsigned patt3[20]
Definition: MdcxSeg.h:9
double Xline_slope()
Definition: MdcxSeg.h:17
double Xline_bbrrf()
Definition: MdcxSeg.h:15
double Phi0_sl_approx()
Definition: MdcxSeg.h:14
double Yline_slope()
Definition: MdcxSeg.h:18
double D0_sl_approx()
Definition: MdcxSeg.h:13
double Yline_bbrrf()
Definition: MdcxSeg.h:16
int Pat()
Definition: MdcxSeg.h:19
virtual ~MdcxTrackFinder()
MdcxTrackFinder(const std::string &name, ISvcLocator *pSvcLocator)
StatusCode finalize()
StatusCode beginRun()
StatusCode execute()
StatusCode initialize()
static void readMCppTable(std::string filenm)
Definition: Pdt.cxx:291
MdcDigiVec & getMdcDigiVec(uint32_t control=0)
static double MdcTime(int timeChannel)
Definition: RawDataUtil.h:8
static double MdcCharge(int chargeChannel)
Definition: RawDataUtil.h:11
const int getFlagLR(void) const
Definition: RecMdcHit.h:47
const double getZhit(void) const
Definition: RecMdcHit.h:55
const int getTrkId(void) const
Definition: RecMdcHit.h:41
const int getId(void) const
Definition: RecMdcHit.h:40
const int getStat(void) const
Definition: RecMdcHit.h:48
const Identifier getMdcId(void) const
Definition: RecMdcHit.h:49
const double getTdc(void) const
Definition: RecMdcHit.h:50
const double getDriftT(void) const
Definition: RecMdcHit.h:52
const double getEntra(void) const
Definition: RecMdcHit.h:54
const double getDriftDistLeft(void) const
Definition: RecMdcHit.h:42
const double getFltLen(void) const
Definition: RecMdcHit.h:56
const double getDoca(void) const
Definition: RecMdcHit.h:53
const HitRefVec getVecHits(void) const
Definition: RecMdcTrack.h:60
const int getNhits() const
Definition: RecMdcTrack.h:49
const double getFiTerm() const
Definition: RecMdcTrack.h:52
virtual void getInfo(double fltLen, HepPoint3D &pos, Hep3Vector &direction) const =0
virtual Hep3Vector momentum(double fltL=0.) const =0
virtual double pt(double fltL=0.) const =0
virtual double chisq() const =0
virtual int charge() const =0
virtual int nDof() const =0
virtual const TrkDifTraj & traj() const =0
virtual HepPoint3D position(double fltL) const =0
void print(std::ostream &ostr) const
Definition: TrkErrCode.cxx:79
int success() const
Definition: TrkErrCode.h:62
double phi0() const
double z0() const
double d0() const
double tanDip() const
virtual void addHistory(const TrkErrCode &status, const char *modulename)
Definition: TrkFit.h:23
virtual int nActive() const =0
virtual TrkExchangePar helix(double fltL) const =0
static bool m_debug
static double nSigmaCut[43]
TrkErrCode fit()
Definition: TrkHitList.cxx:59
hot_iterator begin() const
Definition: TrkHitList.h:45
unsigned nHit() const
Definition: TrkHitList.h:44
TrkHitOnTrk * appendHot(const TrkHitOnTrk *theHot)
Definition: TrkHitList.cxx:91
const TrkHotList & hotList() const
Definition: TrkHitList.cxx:51
hot_iterator end() const
Definition: TrkHitList.h:46
bool removeHit(const TrkFundHit *theHit)
Definition: TrkHitList.cxx:69
void setActivity(bool turnOn)
Definition: TrkHitOnTrk.cxx:96
double fltLen() const
Definition: TrkHitOnTrk.h:91
void setUsability(int usability)
double resid(bool exclude=false) const
double hitRms() const
Definition: TrkHitOnTrk.h:89
const TrkRep * getParentRep() const
Definition: TrkHitOnTrk.h:73
bool isActive() const
Definition: TrkHitOnTrk.h:200
hot_iterator end() const
Definition: TrkHotList.h:45
hot_iterator begin() const
Definition: TrkHotList.h:44
TrkHitList * hits()
Definition: TrkRecoTrk.h:107
virtual void printAll(std::ostream &) const
Definition: TrkRecoTrk.cxx:195
const TrkFit * fitResult() const
Definition: TrkRecoTrk.cxx:387
const TrkFitStatus * status() const
Definition: TrkRecoTrk.cxx:400
virtual double arrivalTime(double fltL) const
Definition: TrkRep.cxx:192
bool setFlipAndDrop(TrkRecoTrk &, bool allowFlips, bool allowDrops) const
TrkRecoTrk * makeTrack(const TrkExchangePar &helix, const double chi2, const TrkContext &, double trackT0) const
_EXTERN_ std::string RecMdcTrackCol
Definition: EventModel.h:86
_EXTERN_ std::string RecMdcHitCol
Definition: EventModel.h:85
char * c_str(Index i)
Definition: EvtCyclic3.cc:252
NTuple::Item< long > g_eventNo
Definition: ntupleItem.h:22
int nhits