00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026 #include "MdcxReco/MdcxTrackFinder.h"
00027
00028
00029
00030
00031 #include <assert.h>
00032 #include <time.h>
00033 #include <math.h>
00034
00035
00036
00037
00038 #include <iostream>
00039
00040
00041
00042
00043 #include "GaudiKernel/MsgStream.h"
00044 #include "GaudiKernel/AlgFactory.h"
00045 #include "BField/BField.h"
00046 #include "MdcRawEvent/MdcDigi.h"
00047 #include "MdcGeom/MdcDetector.h"
00048 #include "EventModel/EventHeader.h"
00049
00050 #include "MdcxReco/MdcxHit.h"
00051 #include "MdcxReco/MdcxHits.h"
00052 #include "MdcxReco/MdcxFindSegs.h"
00053 #include "MdcxReco/MdcxSeg.h"
00054 #include "MdcxReco/MdcxHel.h"
00055 #include "MdcxReco/MdcxFittedHel.h"
00056 #include "MdcxReco/MdcxFindTracks.h"
00057 #include "MdcxReco/Mdcxprobab.h"
00058 #include "CLHEP/Alist/AIterator.h"
00059 #include "TrkBase/TrkExchangePar.h"
00060 #include "TrkBase/TrkErrCode.h"
00061 #include "TrkBase/TrkRecoTrk.h"
00062 #include "TrkBase/TrkFit.h"
00063 #include "TrkBase/TrkHitList.h"
00064 #include "TrkFitter/TrkHelixMaker.h"
00065 #include "TrkFitter/TrkLineMaker.h"
00066 #include "MdcxReco/MdcxAddHits.h"
00067 #include "MdcxReco/MdcxMergeDups.h"
00068 #include "TrkBase/TrkFitStatus.h"
00069 #include "MdcRecEvent/RecMdcTrack.h"
00070 #include "MdcRecEvent/RecMdcHit.h"
00071 #include "ReconEvent/ReconEvent.h"
00072 #include "Identifier/MdcID.h"
00073 #include "Identifier/Identifier.h"
00074 #include "EvTimeEvent/RecEsTime.h"
00075 #include "McTruth/McParticle.h"
00076 #include "McTruth/MdcMcHit.h"
00077 #include "MdcCalibFunSvc/IMdcCalibFunSvc.h"
00078 #include "RawDataProviderSvc/IRawDataProviderSvc.h"
00079 #include "RawDataProviderSvc/MdcRawDataProvider.h"
00080 #include "RawEvent/RawDataUtil.h"
00081 #include "MdcRecoUtil/Pdt.h"
00082 #include "MdcTrkRecon/MdcTrack.h"
00083
00084 #include "TrkBase/TrkHotList.h"
00085 #include "GaudiKernel/INTupleSvc.h"
00086 #include "GaudiKernel/NTuple.h"
00087 #include "MdcxReco/MdcxHistItem.h"
00088 #include "MdcGeom/Constants.h"
00089 #include "MdcGeom/MdcTrkReconCut.h"
00090 #include "MdcxReco/MdcxSegPatterns.h"
00091 #include "TrigEvent/TrigData.h"
00092
00093 #include "BesTimerSvc/IBesTimerSvc.h"
00094 #include "BesTimerSvc/BesTimerSvc.h"
00095 #include "BesTimerSvc/BesTimer.h"
00096
00097
00098
00099 using std::cout;
00100 using std::endl;
00101
00102 extern double MdcTrkReconCut_helix_fit[43];
00103
00104
00105
00106
00107 MdcxTrackFinder::MdcxTrackFinder(const std::string& name, ISvcLocator* pSvcLocator) :
00108 Algorithm(name, pSvcLocator), m_mdcCalibFunSvc(0)
00109 {
00110
00111 declareProperty("pdtFile", m_pdtFile = "pdt.table");
00112
00113 declareProperty("debug", m_debug= 0);
00114 declareProperty("hist", m_hist= 0);
00115 declareProperty("mcHist", m_mcHist = false);
00116
00117 declareProperty("cresol", m_cresol = 0.013);
00118
00119 declareProperty("getDigiFlag", m_getDigiFlag = 0);
00120 declareProperty("maxMdcDigi", m_maxMdcDigi= 0);
00121 declareProperty("keepBadTdc", m_keepBadTdc= 0);
00122 declareProperty("dropHot", m_dropHot= 0);
00123 declareProperty("keepUnmatch", m_keepUnmatch= 0);
00124 declareProperty("salvageTrk", m_salvageTrk = false);
00125 declareProperty("dropMultiHotInLayer",m_dropMultiHotInLayer = false);
00126 declareProperty("dropTrkPt", m_dropTrkPt = -999.);
00127 declareProperty("d0Cut", m_d0Cut = 999.);
00128 declareProperty("z0Cut", m_z0Cut = 999.);
00129
00130 declareProperty("minMdcDigi", m_minMdcDigi = 0);
00131 declareProperty("countPropTime",m_countPropTime = true);
00132 declareProperty("addHitCut", m_addHitCut = 5.);
00133 declareProperty("dropHitsSigma",m_dropHitsSigma);
00134 declareProperty("helixFitCut", m_helixFitCut);
00135 declareProperty("minTrkProb", m_minTrkProb = 0.01);
00136 declareProperty("csmax4", m_csmax4 = 50.);
00137 declareProperty("csmax3", m_csmax3 = 1.);
00138 declareProperty("helixFitSigma",m_helixFitSigma= 5.);
00139 declareProperty("maxRcsInAddSeg",m_maxRcsInAddSeg= 50.);
00140 declareProperty("nSigAddHitTrk", m_nSigAddHitTrk= 5.);
00141 declareProperty("maxProca", m_maxProca= 0.6);
00142 declareProperty("doSag", m_doSag= false);
00143 declareProperty("lineFit", m_lineFit = false);
00144
00145 }
00146
00147
00148
00149
00150 MdcxTrackFinder::~MdcxTrackFinder() {
00151 delete m_bfield;
00152 }
00153
00154 StatusCode MdcxTrackFinder::beginRun(){
00155
00156 m_gm = MdcDetector::instance(m_doSag);
00157 if(NULL == m_gm) return StatusCode::FAILURE;
00158 MdcxHit::setMdcDetector(m_gm);
00159
00160 return StatusCode::SUCCESS;
00161 }
00162
00163
00164
00165 StatusCode MdcxTrackFinder::initialize(){
00166 MsgStream log(msgSvc(), name());
00167 log << MSG::INFO << "in initialize()" << endreq;
00168
00169 t_nTkTot = 0;
00170 for (int i=0;i<20;i++) t_nTkNum[i]=0;
00171
00172
00173 #ifdef MDCXTIMEDEBUG
00174 StatusCode tsc = service( "BesTimerSvc", m_timersvc);
00175 if( tsc.isFailure() ) {
00176 log << MSG::WARNING << name() << ": Unable to locate BesTimer Service" << endreq;
00177 return StatusCode::FAILURE;
00178 }
00179 m_timer[0] = m_timersvc->addItem("Execution");
00180 m_timer[0]->propName("Execution");
00181 #endif
00182
00183 if (m_helixFitCut.size() == 43){
00184 for(int i=0; i<43; i++){
00185
00186 TrkHelixFitter::nSigmaCut[i] = m_helixFitCut[i];
00187 }
00188 }
00189 MdcxParameters::debug = m_debug;
00190 MdcxParameters::minTrkProb = m_minTrkProb;
00191 MdcxParameters::csmax4 = m_csmax4;
00192 MdcxParameters::csmax3 = m_csmax3;
00193 MdcxParameters::helixFitSigma = m_helixFitSigma;
00194 MdcxParameters::maxRcsInAddSeg= m_maxRcsInAddSeg;
00195 MdcxParameters::nSigAddHitTrk = m_nSigAddHitTrk;
00196 MdcxParameters::maxProca = m_maxProca;
00197 TrkHelixFitter::m_debug = (m_debug>7);
00198 Pdt::readMCppTable(m_pdtFile);
00199 MdcxFittedHel::debug = m_debug;
00200
00201
00202
00203 IMdcCalibFunSvc* imdcCalibSvc;
00204 StatusCode sc = service ("MdcCalibFunSvc", imdcCalibSvc);
00205 m_mdcCalibFunSvc = dynamic_cast<MdcCalibFunSvc*>(imdcCalibSvc);
00206 if ( sc.isFailure() ){
00207 log << MSG::FATAL << "Could not load MdcCalibFunSvc!" << endreq;
00208 return StatusCode::FAILURE;
00209 }
00210 MdcxHit::setMdcCalibFunSvc(m_mdcCalibFunSvc);
00211 MdcxHit::setCountPropTime(m_countPropTime);
00212
00213
00214 IRawDataProviderSvc* iRawDataProvider;
00215 sc = service ("RawDataProviderSvc", iRawDataProvider);
00216 if ( sc.isFailure() ){
00217 log << MSG::FATAL << "Could not load RawDataProviderSvc!" << endreq;
00218 return StatusCode::FAILURE;
00219 }
00220 m_rawDataProviderSvc = dynamic_cast<RawDataProviderSvc*>(iRawDataProvider);
00221
00222
00223
00224 sc = service ("MagneticFieldSvc",m_pIMF);
00225 if(sc!=StatusCode::SUCCESS) {
00226 log << MSG::ERROR << "Unable to open Magnetic field service"<<endreq;
00227 }
00228 m_bfield = new BField(m_pIMF);
00229 log << MSG::INFO << "field z = "<<m_bfield->bFieldNominal()<< endreq;
00230 m_context = new TrkContextEv(m_bfield);
00231
00232 if (m_hist) {bookNTuple();}
00233 if (m_dropHitsSigma.size()==43){
00234 for (int ii=0;ii<43;ii++) {
00235 MdcxParameters::dropHitsSigma[ii]=m_dropHitsSigma[ii];
00236 }
00237 }
00238
00239 return StatusCode::SUCCESS;
00240 }
00241
00242 StatusCode MdcxTrackFinder::execute() {
00243
00244
00245 b_saveEvent=false;
00246 setFilterPassed(b_saveEvent);
00247 #ifdef MDCXTIMEDEBUG
00248 m_timer[0]->start();
00249 #endif
00250 MsgStream log(msgSvc(), name());
00251 log << MSG::INFO << "in execute()" << endreq;
00252 StatusCode sc;
00253
00254 nTk = 0; t_nTdsTk = 0;
00255
00256
00257
00258 SmartDataPtr<Event::EventHeader> evtHead(eventSvc(),"/Event/EventHeader");
00259 if (!evtHead) {
00260 log << MSG::FATAL<< "Could not retrieve event header" << endreq;
00261 return StatusCode::FAILURE;
00262 }
00263 m_eventNo = evtHead->eventNumber();
00264
00265 long t_evtNo = m_eventNo;
00266 g_eventNo = m_eventNo;
00267
00268 IDataManagerSvc *dataManSvc;
00269 DataObject *aTrackCol;
00270 DataObject *aHitCol;
00271 if (!m_salvageTrk) {
00272 SmartIF<IDataManagerSvc> dataManSvc(eventSvc());
00273 eventSvc()->findObject("/Event/Recon/RecMdcTrackCol",aTrackCol);
00274 if(aTrackCol != NULL) {
00275 dataManSvc->clearSubTree("/Event/Recon/RecMdcTrackCol");
00276 eventSvc()->unregisterObject("/Event/Recon/RecMdcTrackCol");
00277 }
00278 eventSvc()->findObject("/Event/Recon/RecMdcHitCol",aHitCol);
00279 if(aHitCol != NULL) {
00280 dataManSvc->clearSubTree("/Event/Recon/RecMdcHitCol");
00281 eventSvc()->unregisterObject("/Event/Recon/RecMdcHitCol");
00282 }
00283 }
00284
00285
00286
00287
00288 DataObject *aReconEvent;
00289 eventSvc()->findObject("/Event/Recon",aReconEvent);
00290 if (aReconEvent==NULL) {
00291 aReconEvent = new ReconEvent();
00292 sc = eventSvc()->registerObject("/Event/Recon",aReconEvent);
00293 if(sc != StatusCode::SUCCESS) {
00294 log << MSG::FATAL << "Could not register ReconEvent" <<endreq;
00295 return StatusCode::FAILURE;
00296 }
00297 }
00298 RecMdcTrackCol* trackList;
00299 eventSvc()->findObject("/Event/Recon/RecMdcTrackCol",aTrackCol);
00300 if (aTrackCol) {
00301 trackList = dynamic_cast<RecMdcTrackCol*> (aTrackCol);
00302 }else{
00303 trackList = new RecMdcTrackCol;
00304 sc = eventSvc()->registerObject(EventModel::Recon::RecMdcTrackCol, trackList);
00305 if(!sc.isSuccess()) {
00306 log << MSG::FATAL << " Could not register RecMdcTrack collection" <<endreq;
00307 return StatusCode::FAILURE;
00308 }
00309 }
00310 RecMdcHitCol* hitList;
00311 eventSvc()->findObject("/Event/Recon/RecMdcHitCol",aHitCol);
00312 if (aHitCol) {
00313 hitList = dynamic_cast<RecMdcHitCol*> (aHitCol);
00314 }else{
00315 hitList = new RecMdcHitCol;
00316 sc = eventSvc()->registerObject(EventModel::Recon::RecMdcHitCol, hitList);
00317 if(!sc.isSuccess()) {
00318 log << MSG::FATAL << " Could not register RecMdcHit collection" <<endreq;
00319 return StatusCode::FAILURE;
00320 }
00321 }
00322
00323
00324
00325
00326 DataObject *pnode = 0;
00327 sc = eventSvc()->retrieveObject("/Event/Hit",pnode);
00328 if(!sc.isSuccess()) {
00329 pnode = new DataObject;
00330 sc = eventSvc()->registerObject("/Event/Hit",pnode);
00331 if(!sc.isSuccess()) {
00332 log << MSG::FATAL << " Could not register /Event/Hit branch " <<endreq;
00333 return StatusCode::FAILURE;
00334 }
00335 }
00336 SmartDataPtr<MdcHitCol> m_hitCol(eventSvc(), "/Event/Hit/MdcHitCol");
00337 if (!m_hitCol){
00338 m_hitCol= new MdcHitCol;
00339 sc = eventSvc()->registerObject("/Event/Hit/MdcHitCol",m_hitCol);
00340 if(!sc.isSuccess()) {
00341 log << MSG::FATAL << " Could not register hit collection" <<endreq;
00342 return StatusCode::FAILURE;
00343 }
00344 }
00345
00346
00347
00348
00349
00350 m_bunchT0 = -999.;
00351 SmartDataPtr<RecEsTimeCol> aevtimeCol(eventSvc(),"/Event/Recon/RecEsTimeCol");
00352 if (!aevtimeCol || aevtimeCol->size()==0) {
00353 log << MSG::WARNING<< "evt "<<m_eventNo<<" Could not find RecEsTimeCol"<< endreq;
00354 return StatusCode::SUCCESS;
00355 }
00356
00357 RecEsTimeCol::iterator iter_evt = aevtimeCol->begin();
00358 for(; iter_evt!=aevtimeCol->end(); iter_evt++){
00359 m_bunchT0 = (*iter_evt)->getTest();
00360 m_t0Stat = (*iter_evt)->getStat();
00361 if ((m_t0Stat==0) || (m_bunchT0 < 0.) || (m_bunchT0 > 9999.0) ){
00362 log << MSG::WARNING << "Skip evt:"<<m_eventNo<< " by t0 = "<<m_bunchT0 << endreq;
00363
00364 }
00365 }
00366 if(m_debug>1) std::cout<<name()<<" t0 "<<m_bunchT0<<" t0Stat "<<m_t0Stat<<std::endl;
00367 int trigtiming=-10;
00368 SmartDataPtr<TrigData> trigData(eventSvc(),"/Event/Trig/TrigData");
00369 if(trigData){
00370 log << MSG::INFO <<"Trigger conditions 0--43:"<<endreq;
00371 for(int i = 0; i < 48; i++) {
00372 log << MSG::INFO << trigData->getTrigCondName(i)<<" ---- "<<trigData->getTrigCondition(i)<< endreq;
00373 }
00374 for(int i = 0; i < 16; i++) log << MSG::INFO << "Trigger channel "<< i << ": " << trigData->getTrigChannel(i) << endreq;
00375 m_timing=trigData->getTimingType();
00376
00377 log << MSG::INFO <<"Tigger Timing type: "<< trigtiming << endreq;
00378 }
00379
00380
00381
00382
00383 m_mdcxHits.reset();
00384 uint32_t getDigiFlag = 0;
00385 getDigiFlag += m_maxMdcDigi;
00386 if(m_dropHot||m_salvageTrk) getDigiFlag |= MdcRawDataProvider::b_dropHot;
00387 if(m_keepBadTdc) getDigiFlag |= MdcRawDataProvider::b_keepBadTdc;
00388 if(m_keepUnmatch) getDigiFlag |= MdcRawDataProvider::b_keepUnmatch;
00389 mdcDigiVec = m_rawDataProviderSvc->getMdcDigiVec(getDigiFlag);
00390 t_nDigi = mdcDigiVec.size();
00391
00392 if (0 == t_nDigi ){
00393 log << MSG::WARNING << " No hits in MdcDigiVec" << endreq;
00394 return StatusCode::SUCCESS;
00395 }
00396
00397
00398
00399
00400
00401
00402 if (t_nDigi < m_minMdcDigi){
00403 log << MSG::WARNING << " Skip this event for MdcDigiVec.size() < "<<m_minMdcDigi << endreq;
00404 return StatusCode::SUCCESS;
00405 }
00406 m_mdcxHits.create(mdcDigiVec, m_bunchT0, m_cresol);
00407 const HepAList<MdcxHit>& dchitlist = m_mdcxHits.GetMdcxHitList();
00408
00409
00410
00411
00412 MdcxFindSegs dcsegs(dchitlist,m_debug);
00413 const HepAList<MdcxSeg>& seglist = dcsegs.GetMdcxSeglist();
00414 if(m_debug > 1 ){ dumpMdcxSegs(seglist);}
00415
00416
00417
00418
00419
00420 MdcxFindTracks dctrks(seglist,m_debug);
00421 HepAList<MdcxFittedHel>& firsttrkl = (HepAList<MdcxFittedHel>&)dctrks.GetMdcxTrklist();
00422 if(m_debug>1) dumpTrackList(firsttrkl);
00423
00424
00425
00426
00427
00428
00429
00430
00431
00432
00433 MdcxMergeDups dcmergeem(firsttrkl,m_debug);
00434 HepAList<MdcxFittedHel>& trkl = (HepAList<MdcxFittedHel>&)dcmergeem.GetMergedTrklist();
00435
00436
00437
00438
00439
00440
00441
00442
00443
00444
00445
00446
00447
00448
00449
00450 sc = FitMdcxTrack(trkl, dchitlist, m_hitCol, trackList, hitList);
00451 if (!sc.isSuccess()) {return StatusCode::SUCCESS;}
00452 t_nTdsTk = trackList->size();
00453
00454 t_nTkTot += trackList->size();
00455 if(t_nTdsTk<20) t_nTkNum[t_nTdsTk]++;
00456
00457 #ifdef MDCXTIMEDEBUG
00458 m_timer[0]->stop();
00459 m_xt4time = m_timer[0]->elapsed();
00460 m_xtupleEvt->write();
00461 #endif
00462
00463 if(m_hist) fillEvent();
00464 if (m_debug > 0) {
00465 DataObject* pNode;
00466 eventSvc()->retrieveObject("/Event/Recon/RecMdcTrackCol",pNode);
00467 RecMdcTrackCol *tmpTrackCol = dynamic_cast<RecMdcTrackCol*> (pNode);
00468 eventSvc()->retrieveObject("/Event/Recon/RecMdcHitCol",pNode);
00469 int nTdsTk = 0;
00470 if(tmpTrackCol) nTdsTk = tmpTrackCol->size();
00471
00472
00473 std::cout<< "MdcxTrackFinder: evtNo "<< m_eventNo << " t0="<<m_bunchT0
00474 <<" Found " <<trkl.length()
00475 <<" keep "<< t_nTdsTk
00476 <<" finialy keep "<< nTdsTk;
00477
00478 int ndelete =0; trkl.length() - trackList->size();
00479 if( ndelete>0 ) std::cout <<" delete "<< ndelete;
00480 std::cout <<" track(s)" <<endl;
00481
00482
00483 if(m_debug>1)dumpTdsTrack(tmpTrackCol);
00484 if(m_debug>1)dumpTrack(tmpTrackCol);
00485
00486
00487 }
00488 if((trackList->size()!=4) ) b_saveEvent = true;
00489 setFilterPassed(b_saveEvent);
00490 return StatusCode::SUCCESS;
00491 }
00492
00493 StatusCode MdcxTrackFinder::finalize() {
00494 MsgStream log(msgSvc(), name());
00495 log << MSG::INFO << "in finalize()" << endreq;
00496
00497 std::cout<< " --after " << name() << " keep " << t_nTkTot<< " tracks "<<std::endl;
00498 for(int i=0; i<20; i++){
00499 if(t_nTkNum[i]>0)std::cout<< " nTk="<< i << " "<< t_nTkNum[i]<< std::endl;
00500 }
00501
00502
00503 return StatusCode::SUCCESS;
00504 }
00505
00506
00507 void MdcxTrackFinder::MdcxHitsToHots(MdcxHel& mdcxHelix, const HepAList<MdcxHit>& mdcxHits,
00508 TrkHitList* m_trkHitList, MdcHitCol* mdcHitCol) {
00509
00510 if ( 0 == mdcxHits.length() ) return;
00511
00512 int ihits = 0;
00513 while(mdcxHits[ihits]) {
00514 int ambig = 0;
00515 const MdcHit* newhit = mdcxHits[ihits]->getMdcHit();
00516 if ( 0 == newhit ) {
00517 const MdcDigi* theDigi = mdcxHits[ihits]->getDigi();
00518 int layer = MdcID::layer(mdcxHits[ihits]->getDigi()->identify());
00519 int wire = MdcID::wire(mdcxHits[ihits]->getDigi()->identify());
00520 m_digiMap[layer][wire] = mdcxHits[ihits]->getDigi();
00521 MdcHit *thehit = new MdcHit(theDigi, m_gm);
00522 thehit->setCalibSvc(m_mdcCalibFunSvc);
00523 thehit->setCountPropTime(m_countPropTime);
00524
00525
00526 mdcHitCol->push_back(thehit);
00527 newhit = thehit;
00528 }
00529 MdcRecoHitOnTrack temp(*newhit, ambig, m_bunchT0);
00530 MdcHitOnTrack* newhot = &temp;
00531
00532
00533
00534
00535
00536
00537
00538 m_trkHitList->appendHot(newhot);
00539 newhot->setUsability(true);
00540 newhot->setActivity(true);
00541
00542 ihits++;
00543 }
00544 }
00545
00546 StatusCode MdcxTrackFinder::FitMdcxTrack(HepAList<MdcxFittedHel>& trkl,
00547 const HepAList<MdcxHit>& dchitlist, MdcHitCol* m_hitCol,
00548 RecMdcTrackCol* trackList, RecMdcHitCol* hitList){
00549 StatusCode sc;
00550 MsgStream log(msgSvc(), name());
00551
00552
00553
00554
00555
00556
00557
00558
00559
00560 MdcxAddHits dcaddem(trkl, dchitlist, m_addHitCut);
00561
00562
00563
00564
00565
00566
00567
00568
00569
00570 TrkLineMaker linefactory;
00571 TrkHelixMaker helixfactory;
00572
00573
00574 int tkLen = trkl.length();
00575 for (int kk=0; kk< tkLen; kk++){
00576 const HepAList<MdcxHit>& xhits = trkl[kk]->XHitList();
00577 if(m_debug>2){
00578 std::cout<<__FILE__<<" FitMdcxTrack "<<kk<< std::endl;
00579 for(int i=0; i<xhits.length(); i++){ xhits[i]->print(std::cout); }
00580 std::cout<<std::endl;
00581 }
00582
00583 if(m_debug>2) std::cout<<" before add hits nhits="<<xhits.length()<< std::endl;
00584 HepAList<MdcxHit> xass = dcaddem.GetAssociates(kk);
00585 if(m_debug>2){
00586 std::cout<<" after,add "<<xass.length()<<" hits"<<std::endl;
00587 }
00588
00589 MdcxFittedHel mdcxHelix = *trkl[kk];
00590 double thechisq = mdcxHelix.Chisq();
00591 TrkExchangePar tt(-mdcxHelix.D0(),mdcxHelix.Phi0(),mdcxHelix.Omega(),-mdcxHelix.Z0(),-mdcxHelix.Tanl());
00592 TrkRecoTrk* aTrack;
00593 int nparm;
00594
00595 if (m_lineFit) {
00597 aTrack = linefactory.makeTrack(tt, thechisq, *m_context, m_bunchT0*1.e-9);
00598 nparm = 4;
00599 linefactory.setFlipAndDrop(*aTrack, true, true);
00600 } else {
00601 aTrack = helixfactory.makeTrack(tt, thechisq, *m_context, m_bunchT0*1.e-9);
00602 nparm = 5;
00603 helixfactory.setFlipAndDrop(*aTrack, true, true);
00604 }
00605
00606 TrkHitList* m_trkHitList = aTrack->hits();
00607 if (0 == m_trkHitList) {
00608 delete aTrack;
00609 aTrack = NULL;
00610 continue;
00611 }
00612
00613
00614 MdcxHitsToHots(mdcxHelix, xhits, m_trkHitList, m_hitCol);
00615
00616
00617 MdcxHitsToHots(mdcxHelix, xass, m_trkHitList, m_hitCol);
00618
00619
00620
00621
00622 if(m_dropMultiHotInLayer) dropMultiHotInLayer(m_trkHitList);
00623
00624
00625 if(m_debug>1){ std::cout<< "== put to official fitter "<<std::endl;}
00626 TrkErrCode err = m_trkHitList->fit();
00627
00628 if(m_debug>1){
00629 std::cout<< "== after official fitter "<<std::endl;
00630 aTrack->printAll(std::cout);
00631 }
00632
00633 const TrkFit* theFit = aTrack->fitResult();
00634 float rcs = 10000.0;
00635
00636 if (theFit) {
00637 int ndof = theFit->nActive()-nparm;
00638 if (ndof > 0) rcs = theFit->chisq()/float(ndof);
00639 if (m_debug>1) {
00640 if (4 == nparm) cout << " TrkLineMaker";
00641 else cout << " TrkHelixMaker";
00642 cout << " trkNo. " << kk << " success " << err.success() << " rcs " << rcs
00643 << " chi2 " << theFit->chisq() << " nactive " << theFit->nActive() << endl;
00644 }
00645 }
00646
00647 if ( (1 == err.success()) && (rcs < 20.0) ) {
00648 if(m_debug>1) std::cout<<"== fit success "<<std::endl;
00649 if (4 == nparm) {
00650 linefactory.setFlipAndDrop(*aTrack, false, false);
00651 } else {
00652 helixfactory.setFlipAndDrop(*aTrack, false, false);
00653 }
00654
00655 if (m_debug>1) { cout << "MdcxTrackFinder: accept a track " << endl; }
00656
00657 aTrack->status()->addHistory(err,name().c_str());
00658 store(aTrack, trackList, hitList);
00659 } else if ( (2 == err.success()) && (rcs < 150.0) ) {
00660 if(m_debug > 1) std::cout<<"== fit success = 2, refit now"<<std::endl;
00661 int nrefit = 0;
00662 while (nrefit++ < 5) {
00663 if (m_debug>1) std::cout << "refit time " << nrefit << std::endl;
00664 err = m_trkHitList->fit();
00665 if (err.success() == 1) break;
00666 }
00667 if (err.success() == 1) {
00668 if (4 == nparm) {
00669 linefactory.setFlipAndDrop(*aTrack, false, false);
00670 } else {
00671 helixfactory.setFlipAndDrop(*aTrack, false, false);
00672 }
00673
00674
00675
00676 aTrack->status()->addHistory(err,name().c_str());
00677 store(aTrack, trackList, hitList);
00678 }
00679 } else {
00680 if (m_debug >1) {
00681 std::cout<<"== fit faild "<<std::endl;
00682 err.print(cout);
00683 cout << endl;
00684 }
00685 delete aTrack;
00686 aTrack = NULL;
00687
00688
00689
00690 if(m_debug>1) std::cout<<"== Fit no good; try a better input helix"<<std::endl;
00691 mdcxHelix.Grow(*trkl[kk],xass);
00692 mdcxHelix.VaryRes();
00693 mdcxHelix.SetChiDofBail(1500);
00694 int fail = mdcxHelix.ReFit();
00695 if(m_debug>1)std::cout<<__FILE__<<" refit fail:"<<fail<< std::endl;
00696 if (!mdcxHelix.Fail()) {
00697 const HepAList<MdcxHit>& bxhits = mdcxHelix.XHitList();
00698 thechisq = mdcxHelix.Chisq();
00699 TrkExchangePar tb(mdcxHelix.D0(),mdcxHelix.Phi0(),mdcxHelix.Omega(),mdcxHelix.Z0(),mdcxHelix.Tanl());
00700 TrkRecoTrk* bTrack;
00701 if (4 == nparm){
00702 bTrack = linefactory.makeTrack(tb,thechisq,*m_context,m_bunchT0*1.e-9);
00703 linefactory.setFlipAndDrop(*bTrack, false, false);
00704 }else{
00705 bTrack = helixfactory.makeTrack(tb,thechisq,*m_context,m_bunchT0*1.e-9);
00706 helixfactory.setFlipAndDrop(*bTrack, false, false);
00707 }
00708 TrkHitList* bhits = bTrack->hits();
00709 if (0 == bhits){delete bTrack; bTrack = NULL; continue;}
00710
00711 MdcxHitsToHots(mdcxHelix, bxhits, bhits, m_hitCol);
00712 TrkErrCode berr = bhits->fit();
00713 const TrkFit* bFit = bTrack->fitResult();
00714 rcs=10000.0;
00715 if (bFit) {
00716 int ndof = bFit->nActive() - nparm;
00717 if (ndof > 0) rcs = bFit->chisq()/float(ndof);
00718 if (m_debug >1) {
00719 if (4 == nparm) cout << " TrkLineMaker";
00720 else cout << " TrkHelixMaker";
00721 cout << " success trkNo. " << kk << " status " << berr.success() << " rcs " << rcs
00722 << " chi2 " << bFit->chisq() << " nactive "<< bFit->nActive() << endl;
00723 }
00724 }
00725 if ( ( 1 == berr.success() ) && ( rcs < 50.0 ) ) {
00726
00727 bTrack->status()->addHistory(berr,name().c_str());
00728 if (m_debug>1) {
00729 cout << "MdcxTrackFinder: accept b track and store to TDS" << endl;
00730 bTrack->printAll(cout);
00731 }
00732 store(bTrack, trackList, hitList);
00733 } else {
00734 if (m_debug>1) {
00735 cout<< " fit failed "<<endl;
00736 berr.print(cout);
00737 cout << endl;
00738 }
00739 if (bTrack!=NULL) { delete bTrack; bTrack = NULL; }
00740 }
00741 }else{
00742
00743 }
00744 }
00745 }
00746 if(m_debug >1) dumpTdsTrack(trackList);
00747 return StatusCode::SUCCESS;
00748
00749 }
00750
00751 void MdcxTrackFinder::store(TrkRecoTrk* aTrack, RecMdcTrackCol *trackList,
00752 RecMdcHitCol *hitList) {
00753 assert (aTrack != NULL);
00754 nTk++;
00755 int trackId = trackList->size();
00756 TrkExchangePar helix = aTrack->fitResult()->helix(0.);
00757 if(m_dropTrkPt>0. && (aTrack->fitResult()->pt()<m_dropTrkPt)) {
00758 if(m_debug>1) std::cout<<"MdcxTrackFinder delete track by pt "
00759 <<aTrack->fitResult()->pt()<<"<ptCut "<<m_dropTrkPt << std::endl;
00760 return;
00761 }
00762
00763 if( ( (fabs(helix.d0())>m_d0Cut) ||( fabs(helix.z0())>m_z0Cut) ) ){
00764 if(m_debug>1) std::cout<< name()<<" delete track by d0 "<<helix.d0()<<">d0Cut "<<m_d0Cut
00765 <<" or z0 "<<helix.z0()<<" >z0Cut "<<m_z0Cut << std::endl;
00766 return;
00767 }
00768
00769 if(m_hist) fillTrack(aTrack);
00770 MdcTrack mdcTrack(aTrack);
00771
00772 int tkStat = 1;
00773 int nHitbefore = hitList->size();
00774
00775 mdcTrack.storeTrack(trackId, trackList, hitList, tkStat);
00776 int nHitAfter = hitList->size();
00777 if (nHitAfter - nHitbefore <10 ) setFilterPassed(true);
00778 }
00779
00780 void MdcxTrackFinder::dumpTrack(RecMdcTrackCol* trackList){
00781 RecMdcTrackCol::iterator i_tk = trackList->begin();
00782 for (;i_tk != trackList->end(); i_tk++) {
00783 printTrack(*(i_tk));
00784 }
00785 }
00786
00787 void MdcxTrackFinder::printTrack(RecMdcTrack* tk){
00788
00789 std::cout<< " MdcTrack Id:"<<tk->trackId() <<" q:"<< tk->charge()<< std::endl;
00790 std::cout<< "dr Fi0 Cpa Dz Tanl Chi2 Ndf nSt FiTerm poca" << std::endl;
00791 std::cout<<"(" <<setw(5) << tk->helix(0)<<","<< setw(5) << tk->helix(1)<<"," << setw(5) << tk->helix(2) <<","
00792 << setw(5) << tk->helix(3) << ","<< setw(5) << tk->helix(4) <<")"
00793 << setw(5) << tk->chi2() << setw(4) << tk->ndof()
00794 << setw(4) << tk->getNhits() << setw(4) << tk->nster()
00795 << setw(5) << tk->getFiTerm() <<tk->poca()<<std::endl;
00796 std::cout<< " ErrMat "<<tk->err() << std::endl;
00797
00798 std::cout<< "hitId tkId (l,w) fltLen lr dt ddl tdc "
00799 <<"doca entr z tprop stat " << std::endl;
00800
00801 HitRefVec hl = tk->getVecHits();
00802 HitRefVec::iterator it = hl.begin();
00803 for (;it!=hl.end();++it){
00804 RecMdcHit* h = *it;
00805 int layer = MdcID::layer(h->getMdcId());
00806 double _vprop = (layer<8) ? Constants::vpropInner : Constants::vpropOuter;
00807 const MdcLayer* _layerPtr = m_gm->Layer(layer);
00808 double _zlen = _layerPtr->zLength();
00809 double z = h->getZhit();
00810 double tprop;
00811 if (0 == layer%2){
00812 tprop = (0.5*_zlen + z)/_vprop;
00813 }else{
00814 tprop = (0.5*_zlen - z)/_vprop;
00815 }
00816
00817
00818
00819
00820
00821
00822
00823
00824
00825
00826
00827
00829
00830
00831 std::cout<< setw(3) << h->getId() << setw(2) << h->getTrkId() <<
00832 "(" << MdcID::layer(h->getMdcId()) <<"," << MdcID::wire(h->getMdcId()) <<")"<<
00833 setw(10) << h->getFltLen() <<
00834 setw(2) << h->getFlagLR() <<setw(10) << h->getDriftT() <<
00835 setw(12) << h->getDriftDistLeft() <<
00836 setw(8) << h->getTdc() <<
00837 setw(10) << h->getDoca() <<setw(10) << h->getEntra() <<
00838 setw(10) << h->getZhit() << setw(10) << tprop<<
00839 setw(2)<< h->getStat() << std::endl;
00840 }
00841 }
00842
00843 void MdcxTrackFinder::bookNTuple(){
00844 MsgStream log(msgSvc(), name());
00845 StatusCode sc;
00846 if(!m_hist) return;
00847
00848 g_poison = histoSvc()->book( "poison", "poison",43,0,42,288,0,288 );
00849 g_csmax4 = histoSvc()->book( "csmax4", "csmax4",100,0,100 );
00850 g_csmax3 = histoSvc()->book( "csmax3", "csmax3",50000,0,500000 );
00851 g_omegag = histoSvc()->book( "omegag", "omegag",1000 ,0.,1.);
00852 g_dPhiAU = histoSvc()->book( "dPhiAU", "dPhiAU",1000 ,0.,3.5);
00853 g_dPhiAU_0 = histoSvc()->book( "dPhiAU_0", "dPhiAU_0",1000 ,0.,3.5);
00854 g_dPhiAU_1 = histoSvc()->book( "dPhiAU_1", "dPhiAU_1",1000 ,0.,3.5);
00855 g_dPhiAU_5 = histoSvc()->book( "dPhiAU_5", "dPhiAU_5",1000 ,0.,3.5);
00856 g_dPhiAU_7 = histoSvc()->book( "dPhiAU_7", "dPhiAU_7",1000 ,0.,3.5);
00857 g_dPhiAV = histoSvc()->book( "dPhiAV", "dPhiAV",1000 ,0.,3.5);
00858 g_addSegPhi = histoSvc()->book( "addSegPhi", "addSegPhi",1000 ,0.,3.5);
00859 g_dPhiAV_0 = histoSvc()->book( "dPhiAV_0", "dPhiAV_0",1000 ,0.,3.5);
00860 g_dPhiAV_1 = histoSvc()->book( "dPhiAV_1", "dPhiAV_1",1000 ,0.,3.5);
00861 g_dPhiAV_6 = histoSvc()->book( "dPhiAV_6", "dPhiAV_6",1000 ,0.,3.5);
00862 g_dPhiAV_8 = histoSvc()->book( "dPhiAV_8", "dPhiAV_8",1000 ,0.,3.5);
00863 g_trkllmk = histoSvc()->book( "trkllmk", "trkllmk",43,0.,42);
00864 g_trklgood = histoSvc()->book( "trklgood", "trklgood",43,0.,42);
00865 g_trklcircle = histoSvc()->book( "trklcircle", "trklcircle",43,0.,42);
00866 g_trklhelix= histoSvc()->book( "trklhelix", "trklhelix",43,0.,42);
00867 g_trkldrop1= histoSvc()->book( "trkldrop1", "trkldrop1",43,0.,42);
00868 g_trkldrop2= histoSvc()->book( "trkldrop2", "trkldrop2",43,0.,42);
00869 g_trklappend1= histoSvc()->book( "trklappend1", "trklappend1",43,0.,42);
00870 g_trklappend2= histoSvc()->book( "trklappend2", "trklappend2",43,0.,42);
00871 g_trklappend3= histoSvc()->book( "trklappend3", "trklappend3",43,0.,42);
00872
00873 g_trklfirstProb= histoSvc()->book( "trklfirstProb", "trklfirstProb",200,0.,2.);
00874 g_trkltemp= histoSvc()->book( "trkltemp", "trkltemp",200,-3.5,3.5);
00875
00876 g_trklproca= histoSvc()->book( "trklproca", "trklproca",100,0.,5.);
00877 g_trkle= histoSvc()->book( "trkle", "trkle",200,0.,0.025);
00878 g_trkld= histoSvc()->book( "trkld", "trkld",200,-1.2,1.2);
00879 g_trklel= histoSvc()->book( "trklel", "trklel",200,0.,0.025,43,0,43);
00880 g_trkldl= histoSvc()->book( "trkldl", "trkldl",200,-1.2,1.2,43,0,43);
00881 g_trkldoca= histoSvc()->book( "trkldoca", "trkldoca",200,-1.2,1.2);
00882 g_trkllayer= histoSvc()->book( "trkllayer", "trkllayer",43,0,43);
00883 g_dropHitsSigma= histoSvc()->book( "dropHitsSigma", "dropHitsSigma",43,0,43,100,0,11);
00884 g_addHitCut= histoSvc()->book( "addHitCut", "addHitCut",100,0,10);
00885 g_addHitCut2d= histoSvc()->book( "addHitCut2d", "addHitCut2d",43,0,43,100,0,10);
00886
00887
00888 NTuplePtr nt1(ntupleSvc(), "MdcxReco/rec");
00889 if ( nt1 ) { m_xtuple1 = nt1;}
00890 else {
00891 m_xtuple1 = ntupleSvc()->book ("MdcxReco/rec", CLID_ColumnWiseTuple, "MdcxReco reconsturction results");
00892 if ( m_xtuple1 ) {
00893 sc = m_xtuple1->addItem ("t0", m_xt0);
00894 sc = m_xtuple1->addItem ("timing", m_xtiming);
00895 sc = m_xtuple1->addItem ("t0Stat", m_xt0Stat);
00896 sc = m_xtuple1->addItem ("t0Truth", m_xt0Truth);
00897
00898 sc = m_xtuple1->addItem ("nSlay", m_xnSlay);
00899 sc = m_xtuple1->addItem ("p", m_xp);
00900 sc = m_xtuple1->addItem ("pt", m_xpt);
00901 sc = m_xtuple1->addItem ("pz", m_xpz);
00902 sc = m_xtuple1->addItem ("d0", m_xd0);
00903 sc = m_xtuple1->addItem ("phi0", m_xphi0);
00904 sc = m_xtuple1->addItem ("cpa", m_xcpa);
00905 sc = m_xtuple1->addItem ("z0", m_xz0);
00906 sc = m_xtuple1->addItem ("tanl", m_xtanl);
00907 sc = m_xtuple1->addItem ("q", m_xq);
00908 sc = m_xtuple1->addItem ("pocax", m_xpocax);
00909 sc = m_xtuple1->addItem ("pocay", m_xpocay);
00910 sc = m_xtuple1->addItem ("pocaz", m_xpocaz);
00911
00912 sc = m_xtuple1->addItem ("evtNo", m_xevtNo);
00913 sc = m_xtuple1->addItem ("nSt", m_xnSt);
00914 sc = m_xtuple1->addItem ("nAct", m_xnAct);
00915 sc = m_xtuple1->addItem ("nDof", m_xnDof);
00916 sc = m_xtuple1->addItem ("chi2", m_xchi2);
00917 sc = m_xtuple1->addItem ("tkId", m_xtkId);
00918 sc = m_xtuple1->addItem ("nHit", m_xnHit, 0, 7000);
00919 sc = m_xtuple1->addIndexedItem ("resid", m_xnHit, m_xresid);
00920 sc = m_xtuple1->addIndexedItem ("sigma", m_xnHit, m_xsigma);
00921 sc = m_xtuple1->addIndexedItem ("driftD", m_xnHit, m_xdriftD);
00922 sc = m_xtuple1->addIndexedItem ("driftT", m_xnHit, m_xdriftT);
00923 sc = m_xtuple1->addIndexedItem ("doca", m_xnHit, m_xdoca);
00924 sc = m_xtuple1->addIndexedItem ("entra", m_xnHit, m_xentra);
00925 sc = m_xtuple1->addIndexedItem ("ambig", m_xnHit, m_xambig);
00926 sc = m_xtuple1->addIndexedItem ("fltLen", m_xnHit, m_xfltLen);
00927 sc = m_xtuple1->addIndexedItem ("tof", m_xnHit, m_xtof);
00928 sc = m_xtuple1->addIndexedItem ("act", m_xnHit, m_xact);
00929 sc = m_xtuple1->addIndexedItem ("tdc", m_xnHit, m_xtdc);
00930 sc = m_xtuple1->addIndexedItem ("adc", m_xnHit, m_xadc);
00931 sc = m_xtuple1->addIndexedItem ("layer", m_xnHit, m_xlayer);
00932 sc = m_xtuple1->addIndexedItem ("wire", m_xnHit, m_xwire);
00933 sc = m_xtuple1->addIndexedItem ("x", m_xnHit, m_xx);
00934 sc = m_xtuple1->addIndexedItem ("y", m_xnHit, m_xy);
00935 sc = m_xtuple1->addIndexedItem ("z", m_xnHit, m_xz);
00936 } else {
00937 log << MSG::ERROR << " Cannot book N-tuple: MdcxReco/rec" << endmsg;
00938
00939 }
00940 }
00941 NTuplePtr nt4(ntupleSvc(), "MdcxReco/evt");
00942 if ( nt4 ) { m_xtupleEvt = nt4;}
00943 else {
00944 m_xtupleEvt = ntupleSvc()->book ("MdcxReco/evt", CLID_ColumnWiseTuple, "MdcxReco event data");
00945 if ( m_xtupleEvt ) {
00946 sc = m_xtupleEvt->addItem ("evtNo", m_xt4EvtNo );
00947 sc = m_xtupleEvt->addItem ("nRecTk", m_xt4nRecTk );
00948 sc = m_xtupleEvt->addItem ("nTdsTk", m_xt4nTdsTk );
00949 sc = m_xtupleEvt->addItem ("t0", m_xt4t0);
00950 sc = m_xtupleEvt->addItem ("t0Stat", m_xt4t0Stat);
00951 sc = m_xtupleEvt->addItem ("t0Truth", m_xt4t0Truth);
00952 sc = m_xtupleEvt->addItem ("time", m_xt4time);
00953 sc = m_xtupleEvt->addItem ("nDigi", m_xt4nDigi, 0, 7000);
00954 sc = m_xtupleEvt->addIndexedItem ("layer", m_xt4nDigi, m_xt4Layer);
00955 sc = m_xtupleEvt->addIndexedItem ("rt", m_xt4nDigi, m_xt4Time);
00956 sc = m_xtupleEvt->addIndexedItem ("rc", m_xt4nDigi, m_xt4Charge);
00957
00958
00959 } else {
00960 log << MSG::ERROR << "Cannot book N-tuple: MdcxReco/evt" << endmsg;
00961
00962 }
00963 }
00964
00965
00966 NTuplePtr ntSeg(ntupleSvc(), "MdcxReco/seg");
00967 if ( ntSeg ) { m_xtupleSeg = ntSeg;}
00968 else {
00969 m_xtupleSeg = ntupleSvc()->book ("MdcxReco/seg", CLID_ColumnWiseTuple, "MdcxTrackFinder segment data");
00970 if ( m_xtupleSeg ) {
00971 sc = m_xtupleSeg->addItem ("sl", m_xtsSl);
00972 sc = m_xtupleSeg->addItem ("d0", m_xtsD0);
00973 sc = m_xtupleSeg->addItem ("omega", m_xtsOmega);
00974 sc = m_xtupleSeg->addItem ("phi0", m_xtsPhi0);
00975 sc = m_xtupleSeg->addItem ("d0sl", m_xtsD0_sl_approx);
00976 sc = m_xtupleSeg->addItem ("phi0sl", m_xtsPhi0_sl_approx);
00977 sc = m_xtupleSeg->addItem ("xbbrrf", m_xtsXline_bbrrf);
00978 sc = m_xtupleSeg->addItem ("ybbrrf", m_xtsYline_bbrrf);
00979 sc = m_xtupleSeg->addItem ("xslope", m_xtsXline_slope);
00980 sc = m_xtupleSeg->addItem ("yslope", m_xtsYline_slope);
00981 sc = m_xtupleSeg->addItem ("pat", m_xtsPat);
00982 sc = m_xtupleSeg->addItem ("chisq", m_xtsChisq);
00983 sc = m_xtupleSeg->addItem ("nDigi", m_xtsNDigi, 0, 20);
00984 sc = m_xtupleSeg->addIndexedItem ("layer", m_xtsNDigi, m_xtsLayer);
00985 sc = m_xtupleSeg->addIndexedItem ("wire", m_xtsNDigi, m_xtsWire);
00986 sc = m_xtupleSeg->addIndexedItem ("inSeg", m_xtsNDigi, m_xtsInSeg);
00987 }else{
00988 log << MSG::ERROR << "Cannot book N-tuple: MdcxReco/evt" << endmsg;
00989
00990 }
00991 }
00992 NTuplePtr nt5(ntupleSvc(), "MdcxReco/trkl");
00993 if ( nt5 ) { m_xtupleTrkl= nt5;}
00994 else {
00995 m_xtupleTrkl= ntupleSvc()->book ("MdcxReco/trkl", CLID_RowWiseTuple, "MdcxReco track info");
00996 if ( m_xtupleTrkl) {
00997 sc = m_xtupleTrkl->addItem ("layer", m_xt5Layer);
00998 sc = m_xtupleTrkl->addItem ("wire", m_xt5Wire);
00999 }else{
01000 log << MSG::ERROR << "Cannot book N-tuple: MdcxReco/trkl" << endmsg;
01001
01002 }
01003 }
01004
01005 NTuplePtr ntCsmcSew(ntupleSvc(), "MdcxReco/csmc");
01006 if ( ntCsmcSew ) { m_xtupleCsmcSew = ntCsmcSew;}
01007 else {
01008 m_xtupleCsmcSew = ntupleSvc()->book ("MdcxReco/csmc", CLID_ColumnWiseTuple, "MdcxReco reconsturction results");
01009 if ( m_xtupleCsmcSew ) {
01010 sc = m_xtupleCsmcSew->addItem ("dD0", m_csmcD0);
01011 sc = m_xtupleCsmcSew->addItem ("dPhi0", m_csmcPhi0);
01012 sc = m_xtupleCsmcSew->addItem ("dZ0", m_csmcZ0);
01013 sc = m_xtupleCsmcSew->addItem ("dOmega", m_csmcOmega);
01014 sc = m_xtupleCsmcSew->addItem ("dPt", m_csmcPt);
01015 sc = m_xtupleCsmcSew->addItem ("dTanl", m_csmcTanl);
01016 }else{
01017 log << MSG::ERROR << "Cannot book N-tuple: MdcxReco/csmc" << endmsg;
01018
01019 }
01020 }
01021 NTuplePtr ntSegAdd1(ntupleSvc(), "MdcxReco/addSeg1");
01022 if ( ntSegAdd1 ) { m_xtupleAddSeg1 = ntSegAdd1;}
01023 else {
01024 m_xtupleAddSeg1 = ntupleSvc()->book ("MdcxReco/addSeg1", CLID_ColumnWiseTuple, "MdcxReco event data");
01025 if ( m_xtupleAddSeg1 ) {
01026 sc = m_xtupleAddSeg1->addItem ("same", m_addSegSame);
01027 sc = m_xtupleAddSeg1->addItem ("seedSl", m_addSegSeedSl );
01028 sc = m_xtupleAddSeg1->addItem ("seedPhi", m_addSegSeedPhi );
01029 sc = m_xtupleAddSeg1->addItem ("seedPhiLay",m_addSegSeedPhiLay );
01030 sc = m_xtupleAddSeg1->addItem ("seedLen", m_addSegSeedLen );
01031 sc = m_xtupleAddSeg1->addItem ("seedD0", m_addSegSeedD0 );
01032 sc = m_xtupleAddSeg1->addItem ("seedPhi0", m_addSegSeedPhi0 );
01033 sc = m_xtupleAddSeg1->addItem ("addSl", m_addSegAddSl );
01034 sc = m_xtupleAddSeg1->addItem ("addPhi", m_addSegAddPhi );
01035 sc = m_xtupleAddSeg1->addItem ("addPhiLay", m_addSegAddPhiLay );
01036 sc = m_xtupleAddSeg1->addItem ("addLen", m_addSegAddLen );
01037 sc = m_xtupleAddSeg1->addItem ("addD0", m_addSegAddD0 );
01038 sc = m_xtupleAddSeg1->addItem ("addPhi0", m_addSegAddPhi0 );
01039 }
01040 else{
01041 log << MSG::ERROR << "Cannot book N-tuple: MdcxReco/addSeg1" << endmsg;
01042
01043 }
01044 }
01045
01046 NTuplePtr ntSegComb(ntupleSvc(), "MdcxReco/segComb");
01047 if ( ntSegComb ) { m_xtupleSegComb = ntSegComb;}
01048 else {
01049 m_xtupleSegComb = ntupleSvc()->book ("MdcxReco/segComb", CLID_ColumnWiseTuple, "MdcxReco event data");
01050 if ( m_xtupleSegComb ) {
01051 sc = m_xtupleSegComb->addItem ("evtNo", m_segCombEvtNo);
01052 sc = m_xtupleSegComb->addItem ("omega", m_segCombOmega);
01053 sc = m_xtupleSegComb->addItem ("sameAU", m_segCombSameAU);
01054 sc = m_xtupleSegComb->addItem ("sameUV", m_segCombSameUV);
01055 sc = m_xtupleSegComb->addItem ("dlenAU", m_segCombDLenAU);
01056 sc = m_xtupleSegComb->addItem ("dlenUV", m_segCombDLenUV);
01057 sc = m_xtupleSegComb->addItem ("slA", m_segCombSlA);
01058 sc = m_xtupleSegComb->addItem ("slU", m_segCombSlU);
01059 sc = m_xtupleSegComb->addItem ("slV", m_segCombSlV);
01060 sc = m_xtupleSegComb->addItem ("phiA", m_segCombPhiA);
01061 sc = m_xtupleSegComb->addItem ("phiU", m_segCombPhiU);
01062 sc = m_xtupleSegComb->addItem ("phiV", m_segCombPhiV);
01063 }
01064 else{
01065 log << MSG::ERROR << "Cannot book N-tuple: MdcxReco/segComb" << endmsg;
01066
01067 }
01068 }
01069
01070 NTuplePtr ntDropHits(ntupleSvc(), "MdcxReco/dropHits");
01071 if ( ntDropHits ) { m_xtupleDropHits = ntDropHits;}
01072 else {
01073 m_xtupleDropHits = ntupleSvc()->book ("MdcxReco/dropHits", CLID_ColumnWiseTuple, "MdcxReco event data");
01074 if ( m_xtupleDropHits ) {
01075 sc = m_xtupleDropHits->addItem ("evtNo", m_segDropHitsEvtNo);
01076 sc = m_xtupleDropHits->addItem ("layer", m_segDropHitsLayer);
01077 sc = m_xtupleDropHits->addItem ("wire", m_segDropHitsWire);
01078 sc = m_xtupleDropHits->addItem ("pull", m_segDropHitsPull);
01079 sc = m_xtupleDropHits->addItem ("doca", m_segDropHitsDoca);
01080 sc = m_xtupleDropHits->addItem ("sigma", m_segDropHitsSigma);
01081 sc = m_xtupleDropHits->addItem ("drift", m_segDropHitsDrift);
01082 sc = m_xtupleDropHits->addItem ("mcTkId", m_segDropHitsMcTkId);
01083 } else{
01084 log << MSG::ERROR << "Cannot book N-tuple: MdcxReco/dropHits" << endmsg;
01085
01086 }
01087 }
01088 NTuplePtr ntAddSeg2(ntupleSvc(), "MdcxReco/addSeg2");
01089 if ( ntAddSeg2 ) { m_xtupleAddSeg2 = ntAddSeg2;}
01090 else {
01091 m_xtupleAddSeg2 = ntupleSvc()->book ("MdcxReco/addSeg2", CLID_ColumnWiseTuple, "MdcxReco event data");
01092 if ( m_xtupleAddSeg2 ) {
01093 sc = m_xtupleAddSeg2->addItem ("evtNo", m_addSegEvtNo);
01094 sc = m_xtupleAddSeg2->addItem ("slayer", m_addSegSlayer);
01095 sc = m_xtupleAddSeg2->addItem ("poca", m_addSegPoca);
01096 sc = m_xtupleAddSeg2->addItem ("len", m_addSegLen);
01097 } else{
01098 log << MSG::ERROR << "Cannot book N-tuple: MdcxReco/addSeg2" << endmsg;
01099
01100 }
01101 }
01102
01103 }
01104
01105
01106 void MdcxTrackFinder::fillEvent(){
01107
01108 if (t_nDigi<=0) return;
01109 if (m_xtupleEvt == NULL ) return;
01110 m_xt4EvtNo = m_eventNo;
01111 m_xt4t0 = m_bunchT0;
01112 m_xt4t0Stat = m_t0Stat;
01113 m_xt4t0Truth = m_t0Truth;
01114 m_xt4nRecTk = nTk;
01115 m_xt4nTdsTk = t_nTdsTk;
01116 m_xt4nDigi = t_nDigi;
01117 int iDigi=0;
01118 MdcDigiCol::iterator iter = mdcDigiVec.begin();
01119 for (;iDigi<t_nDigi; iter++ ) {
01120 int l = MdcID::layer((*iter)->identify());
01121 m_xt4Layer[iDigi] = l;
01122
01123 m_xt4Time[iDigi] = RawDataUtil::MdcTime((*iter)->getTimeChannel());
01124 m_xt4Charge[iDigi] = RawDataUtil::MdcCharge((*iter)->getChargeChannel());
01125
01126 iDigi++;
01127 }
01128 m_xtupleEvt->write();
01129 }
01130
01131 void MdcxTrackFinder::fillTrack(TrkRecoTrk* atrack){
01132
01133 if(!m_xtuple1) return;
01134
01135 m_xevtNo = m_eventNo;
01136 int recHitMap[43][288]={0};
01137
01138
01139 if (atrack==NULL) return;
01140
01141 const TrkFit* fitResult = atrack->fitResult();
01142 if (fitResult == 0) return;
01143
01144
01145 int nSt=0;
01146 int seg[11] = {0}; int segme;
01147
01148 HitRefVec hitRefVec;
01149 const TrkHitList* hitlist = atrack->hits();
01150
01151 TrkHitList::hot_iterator hot = hitlist->begin();
01152 int layerCount[43]={0};
01153 int i=0;
01154 for (;hot!= hitlist->end();hot++){
01155
01156 const MdcRecoHitOnTrack* recoHot;
01157 recoHot = dynamic_cast<const MdcRecoHitOnTrack*> (&(*hot));
01158 int layer = recoHot->mdcHit()->layernumber();
01159 int wire = recoHot->mdcHit()->wirenumber();
01160 m_xlayer[i] = layer;
01161
01162 m_xwire[i] = wire;
01163 layerCount[layer]++;
01164 recHitMap[layer][wire]++;
01165 m_xambig[i] = recoHot->wireAmbig();
01166
01167 double fltLen = recoHot->fltLen();
01168 m_xfltLen[i] = fltLen;
01169 double tof = recoHot->getParentRep()->arrivalTime(fltLen);
01170
01171 HepPoint3D pos; Hep3Vector dir;
01172 recoHot->getParentRep()->traj().getInfo(fltLen,pos,dir);
01173 m_xx[i] = pos.x();
01174 m_xy[i] = pos.y();
01175 m_xz[i] = pos.z();
01176 m_xdriftT[i] = recoHot->mdcHit()->driftTime(tof,pos.z());
01177 m_xtof[i] = tof;
01178 m_xdriftD[i] = recoHot->drift();
01179 m_xsigma[i] = recoHot->hitRms();
01180 m_xtdc[i] = recoHot->rawTime();
01181 m_xadc[i] = recoHot->mdcHit()->charge();
01182 m_xdoca[i] = recoHot->dcaToWire();
01183 m_xentra[i] = recoHot->entranceAngle();
01184
01185 m_xact[i] = recoHot->isActive();
01186
01187 double res=999.,rese=999.;
01188 if (recoHot->resid(res,rese,false)){
01189 }else{}
01190 m_xresid[i] = res;
01191
01192 segme=0;
01193 if ( layer >0 ) {segme=(layer-1)/4;}
01194 seg[segme]++;
01195 if (recoHot->layer()->view()) { ++nSt; }
01196 i++;
01197 }
01198
01199 int nSlay=0;
01200 for(int i=0;i<11;i++){
01201 if (seg[i]>0) nSlay++;
01202 }
01203 m_xnSlay = nSlay;
01204
01205
01206
01207
01208 double fltLenPoca = 0.0;
01209 TrkExchangePar helix = fitResult->helix(fltLenPoca);
01210 m_xphi0 = helix.phi0();
01211 m_xtanl = helix.tanDip();
01212 m_xz0 = helix.z0();
01213 m_xd0 = helix.d0();
01214 if(fabs(m_xz0)>20||fabs(m_xd0)>2) {
01215
01216 if(m_debug>1) std::cout<<"evt:"<<m_eventNo<<" BigVtx:"
01217 <<" d0 "<<helix.d0()<<" z0 "<<helix.z0()<<std::endl;
01218 }
01219 m_xpt = fitResult->pt();
01220 if (m_xpt > 0.00001){
01221 m_xcpa = fitResult->charge()/fitResult->pt();
01222 }
01223
01224 CLHEP::Hep3Vector p1 = fitResult->momentum(fltLenPoca);
01225 double px,py,pz,pxy;
01226 pxy = fitResult->pt();
01227 px = p1.x();
01228 py = p1.y();
01229 pz = p1.z();
01230 m_xp = p1.mag();
01231 m_xpz = pz;
01232 HepPoint3D poca = fitResult->position(fltLenPoca);
01233 m_xpocax = poca.x();
01234 m_xpocay = poca.y();
01235 m_xpocaz = poca.z();
01236
01237 m_xq = fitResult->charge();
01238 m_xnAct = fitResult->nActive();
01239 m_xnHit = hitlist->nHit();
01240 m_xnDof = fitResult->nDof();
01241 m_xnSt = nSt;
01242 m_xchi2 = fitResult->chisq();
01243
01244 m_xt0 = m_bunchT0;
01245 m_xtiming = m_timing;
01246 m_xt0Stat = m_t0Stat;
01247 m_xt0Truth= m_t0Truth;
01248 m_xtuple1->write();
01249
01250
01251 }
01252
01253 void MdcxTrackFinder::dumpMdcxSegs(const HepAList<MdcxSeg>& segList) const {
01254 cout << name()<<" found " <<segList.length() << " segs :" << endl;
01255 for (int i =0; i< segList.length(); i++){
01256 std::cout<<i<<" ";
01257 segList[i]->printSeg();
01258 }
01259 }
01260
01261 void MdcxTrackFinder::fillMdcxSegs(const HepAList<MdcxSeg>& segList) const {
01262 if(!m_xtupleSeg) return;
01263
01264 int cellMax[43] ={
01265 40,44,48,56, 64,72,80,80, 76,76,88,88,
01266 100,100,112,112, 128,128,140,140, 160,160,160,160,
01267 176,176,176,176, 208,208,208,208, 240,240,240,240,
01268 256,256,256,256, 288,288,288 };
01269
01270 int hitInSegList[43][288];
01271 for (int ii=0;ii<43;ii++){
01272 for (int jj=0;jj<cellMax[ii];jj++){ hitInSegList[ii][jj] = 0; }
01273 }
01274 for (int i = 0; i < segList.length(); i++) {
01275 MdcxSeg* aSeg = segList[i];
01276 for (int ihit=0 ; ihit< aSeg->Nhits() ; ihit++){
01277 const MdcxHit* hit = aSeg->XHitList()[ihit];
01278 int layer = hit->Layer();
01279 int wire = hit->WireNo();
01280 hitInSegList[layer][wire]++;
01281 }
01282 }
01283 for (int i = 0; i < segList.length(); i++) {
01284 MdcxSeg* aSeg = segList[i];
01285 if(aSeg==NULL)continue;
01286 m_xtsSl = aSeg->SuperLayer();
01287 m_xtsD0 = aSeg->D0();
01288 m_xtsOmega = aSeg->Omega();
01289 m_xtsPhi0 = aSeg->Phi0();
01290 m_xtsD0_sl_approx = aSeg->D0_sl_approx();
01291 m_xtsPhi0_sl_approx = aSeg->Phi0_sl_approx();
01292 m_xtsXline_bbrrf = aSeg->Xline_bbrrf();
01293 m_xtsYline_bbrrf = aSeg->Yline_bbrrf();
01294 m_xtsXline_slope = aSeg->Xline_slope();
01295 m_xtsYline_slope = aSeg->Yline_slope();
01296 int patIndex = -1;
01297 for (int ii = 0;ii<14;ii++){
01298 if (aSeg->Pat() == (int) MdcxSegPatterns::patt4[ii]){
01299 patIndex = ii;
01300 break;
01301 }
01302 }
01303 for (int ii = 0;ii<20;ii++){
01304 if (aSeg->Pat() == (int) MdcxSegPatterns::patt3[ii]){
01305 patIndex = ii +15;
01306 break;
01307 }
01308 }
01309 m_xtsPat = patIndex;
01310 m_xtsChisq = aSeg->Chisq();
01311 m_xtsNDigi = aSeg->Nhits();
01312 for (int ihit=0 ; ihit< aSeg->Nhits() ; ihit++){
01313 const MdcxHit* hit = aSeg->XHitList()[ihit];
01314 int layer = hit->Layer();
01315 int wire = hit->WireNo();
01316 m_xtsLayer[ihit] = layer;
01317 m_xtsWire[ihit] = wire;
01318 m_xtsInSeg[ihit] = hitInSegList[layer][wire];
01319
01320 }
01321 m_xtupleSeg->write();
01322 }
01323
01324 }
01325
01326 void MdcxTrackFinder::dumpTrackList(const HepAList<MdcxFittedHel>& trackList) const {
01327 std::cout<< "dump track list " <<std::endl;
01328 for (int i=0;i<trackList.length();i++){
01329 std::cout<< "track "<<i<<std::endl;
01330 for (int ihit=0 ; ihit< trackList[i]->Nhits() ; ihit++){
01331 const MdcxHit* hit = trackList[i]->XHitList()[ihit];
01332 int layer = hit->Layer();
01333 int wire = hit->WireNo();
01334 std::cout<< " ("<<layer << ","<< wire<<") ";
01335 }
01336 std::cout<<std::endl;
01337 }
01338 }
01339
01340 void MdcxTrackFinder::fillTrkl(const HepAList<MdcxFittedHel>& firsttrkl) const {
01341 if(!m_xtupleTrkl) return;
01342 int nDigi = 0;
01343 int iDigi = 0;
01344 for (int i =0; i< firsttrkl.length(); i++){
01345 nDigi+= firsttrkl[i]->Nhits();
01346 }
01347 for (int i=0;i<firsttrkl.length();i++){
01348 for (int ihit=0 ; ihit< firsttrkl[i]->Nhits() ; ihit++){
01349 const MdcxHit* hit = firsttrkl[i]->XHitList()[ihit];
01350 int layer = hit->Layer();
01351 int wire = hit->WireNo();
01352 m_xt5Layer = layer;
01353 m_xt5Wire = wire;
01354 m_xtupleTrkl->write();
01355 iDigi++;
01356 }
01357 }
01358 }
01359
01360 void MdcxTrackFinder::fillMcTruth(){
01361 MsgStream log(msgSvc(), name());
01362 StatusCode sc;
01363
01364 for (int ii=0;ii<43;ii++){
01365 for (int jj=0;jj<288;jj++){
01366 haveDigi[ii][jj] = -2;
01367 }
01368 }
01369 int nDigi = mdcDigiVec.size();
01370 for (int iDigi =0 ;iDigi<nDigi; iDigi++ ) {
01371 int l = MdcID::layer((mdcDigiVec[iDigi])->identify());
01372 int w = MdcID::wire((mdcDigiVec[iDigi])->identify());
01373
01374 haveDigi[l][w]=1;
01375 }
01376
01377 if( m_mcHist ){
01378
01379 m_t0Truth = -10.;
01380 SmartDataPtr<Event::McParticleCol> mcParticleCol(eventSvc(),"/Event/MC/McParticleCol");
01381 if(!mcParticleCol){
01382 log << MSG::INFO << "Could not retrieve McParticelCol" << endreq;
01383 }else {
01384 Event::McParticleCol::iterator iter_mc = mcParticleCol->begin();
01385 for (;iter_mc != mcParticleCol->end(); iter_mc++){
01386 if ((*iter_mc)->primaryParticle()){
01387 m_t0Truth = (*iter_mc)->initialPosition().t();
01388
01389
01390
01391 }
01392 }
01393 }
01394 }
01395
01396
01397
01398
01399
01400
01401
01402
01403
01404
01405
01406
01407
01408
01409
01410
01411
01412
01413
01414
01415 }
01416
01417 void MdcxTrackFinder::dropMultiHotInLayer(TrkHitList* trkHitList){
01418 double tdr[43];
01419 double tdr_wire[43];
01420 for(int i=0; i<43; i++){tdr[i]=9999.;}
01421
01422
01423 TrkHotList::hot_iterator hotIter = trkHitList->hotList().begin();
01424 while (hotIter!=trkHitList->hotList().end()) {
01425 MdcHitOnTrack* hot = const_cast<MdcHitOnTrack*> (&(*hotIter->mdcHitOnTrack()));
01426
01427
01428 double dt = hot->mdcHit()->driftTime(0.,0.);
01429 int layer = hot->mdcHit()->layernumber();
01430 int wire = hot->mdcHit()->wirenumber();
01431
01432
01433 if (dt < tdr[layer]) {
01434 tdr[layer] = dt;
01435 tdr_wire[layer] = wire;
01436 }
01437 hotIter++;
01438 }
01439
01440
01441
01442
01443
01444
01445 hotIter = trkHitList->hotList().begin();
01446 while (hotIter!=trkHitList->hotList().end()) {
01447 int layer = hotIter->mdcHitOnTrack()->mdcHit()->layernumber();
01448 int wire = hotIter->mdcHitOnTrack()->mdcHit()->wirenumber();
01449
01450
01451 if ((tdr[layer] <9998.) && (tdr_wire[layer]!=wire)){
01452 MdcHitOnTrack* hot = const_cast<MdcHitOnTrack*> (&(*hotIter->mdcHitOnTrack()));
01453 hot->setActivity(false);
01454 trkHitList->removeHit( hotIter->mdcHitOnTrack()->mdcHit() );
01455
01456 }else{
01457 hotIter++;
01458 }
01459 }
01460 }
01461 void MdcxTrackFinder::dumpTdsTrack(RecMdcTrackCol* trackList){
01462 std::cout<< "tksize = "<<trackList->size() << std::endl;
01463 RecMdcTrackCol::iterator it = trackList->begin();
01464 for (;it!= trackList->end();it++){
01465 RecMdcTrack *tk = *it;
01466 std::cout<< "//====RecMdcTrack "<<tk->trackId()<<"====:" << std::endl;
01467 cout <<" d0 "<<tk->helix(0)
01468 <<" phi0 "<<tk->helix(1)
01469 <<" cpa "<<tk->helix(2)
01470 <<" z0 "<<tk->helix(3)
01471 <<" tanl "<<tk->helix(4)
01472 <<endl;
01473 std::cout<<" q "<<tk->charge()
01474 <<" theta "<<tk->theta()
01475 <<" phi "<<tk->phi()
01476 <<" x0 "<<tk->x()
01477 <<" y0 "<<tk->y()
01478 <<" z0 "<<tk->z()
01479 <<" r0 "<<tk->r()
01480 <<endl;
01481 std::cout <<" p "<<tk->p()
01482 <<" pt "<<tk->pxy()
01483 <<" px "<<tk->px()
01484 <<" py "<<tk->py()
01485 <<" pz "<<tk->pz()
01486 <<endl;
01487 std::cout<<" tkStat "<<tk->stat()
01488 <<" chi2 "<<tk->chi2()
01489 <<" ndof "<<tk->ndof()
01490 <<" nhit "<<tk->getNhits()
01491 <<" nst "<<tk->nster()
01492 <<endl;
01493
01494 int nhits = tk->getVecHits().size();
01495 std::cout<<nhits <<" Hits: " << std::endl;
01496 for(int ii=0; ii <nhits ; ii++){
01497 Identifier id(tk->getVecHits()[ii]->getMdcId());
01498 int layer = MdcID::layer(id);
01499 int wire = MdcID::wire(id);
01500 cout<<"("<< layer <<","<<wire<<","<<tk->getVecHits()[ii]->getStat()
01501 <<",lr:"<<tk->getVecHits()[ii]->getFlagLR()<<") ";
01502 }
01503 std::cout << " "<< std::endl;
01504 }
01505 std::cout << " "<< std::endl;
01506
01507
01508 }
01509
01510 void MdcxTrackFinder::dumpTdsHits(RecMdcHitCol* hitList){
01511 std::cout<<__FILE__<<" "<<__LINE__<<" All hits in TDS, nhit="<<hitList->size()<< std::endl;
01512 RecMdcHitCol::iterator it = hitList->begin();
01513 for(;it!= hitList->end(); it++){
01514 RecMdcHit* h = (*it);
01515 Identifier id(h->getMdcId());
01516 int layer = MdcID::layer(id);
01517 int wire = MdcID::wire(id);
01518 cout<<"("<< layer <<","<<wire<<") lr:"<<h->getFlagLR()<<" stat:"<<h->getStat()<<" tk:"<<h->getTrkId() <<" doca:"<<setw(10)<<h->getDoca()<<std::endl;
01519 }
01520 }