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/MdcxCosmicSewer.h"
00027
00028
00029
00030
00031 #include "GaudiKernel/MsgStream.h"
00032 #include "GaudiKernel/AlgFactory.h"
00033 #include "EventModel/EventHeader.h"
00034 #include "EvTimeEvent/RecEsTime.h"
00035 #include "Identifier/MdcID.h"
00036 #include "Identifier/Identifier.h"
00037 #include "GaudiKernel/SmartDataPtr.h"
00038 #include "MdcCalibFunSvc/IMdcCalibFunSvc.h"
00039 #include "RawDataProviderSvc/IRawDataProviderSvc.h"
00040 #include "RawDataProviderSvc/MdcRawDataProvider.h"
00041 #include "MdcRawEvent/MdcDigi.h"
00042 #include "TrkBase/TrkFit.h"
00043 #include "TrkFitter/TrkHelixMaker.h"
00044 #include "TrkFitter/TrkLineMaker.h"
00045 #include "MdcxReco/Mdcxprobab.h"
00046 #include "MdcData/MdcHit.h"
00047 #include "MdcData/MdcRecoHitOnTrack.h"
00048 #include "MdcTrkRecon/MdcTrack.h"
00049 #include "TrkBase/TrkHitList.h"
00050 #include "MdcPrintSvc/IMdcPrintSvc.h"
00051 #include "TrkBase/TrkExchangePar.h"
00052 #include "MdcGeom/EntranceAngle.h"
00053
00054
00055
00056
00057
00058 MdcxCosmicSewer::MdcxCosmicSewer(const std::string& name, ISvcLocator* pSvcLocator) :
00059 Algorithm(name, pSvcLocator)
00060 {
00061
00062 declareProperty("debug", m_debug = false);
00063 declareProperty("hist", m_hist = false);
00064
00065 declareProperty("doSag", m_doSag = false);
00066
00067 declareProperty("cosmicSewPar", m_cosmicSewPar);
00068 declareProperty("cosmicSewSkip", m_cosmicSewSkip=false);
00069
00070 declareProperty("countPropTime", m_countPropTime = true);
00071 declareProperty("lineFit", m_lineFit = false);
00072 }
00073
00074
00075
00076
00077 MdcxCosmicSewer::~MdcxCosmicSewer() {
00078 delete m_bfield;
00079 delete m_context;
00080 }
00081
00082 StatusCode MdcxCosmicSewer::beginRun(){
00083
00084 m_gm = MdcDetector::instance(m_doSag);
00085 if(NULL == m_gm) return StatusCode::FAILURE;
00086
00087 return StatusCode::SUCCESS;
00088 }
00089
00090
00091
00092
00093 StatusCode MdcxCosmicSewer::initialize(){
00094 MsgStream log(msgSvc(), name());
00095 StatusCode sc;
00096 log << MSG::INFO << "in initialize()" << endreq;
00097
00098 i_evt=0;
00099 m_nSewed = 0;
00100
00101
00102 sc = service ("MagneticFieldSvc",m_pIMF);
00103 if(sc!=StatusCode::SUCCESS) {
00104 log << MSG::ERROR << "Unable to open Magnetic field service"<<endreq;
00105 }
00106 m_bfield = new BField(m_pIMF);
00107 log << MSG::INFO << "field z = "<<m_bfield->bFieldNominal()<< endreq;
00108 m_context = new TrkContextEv(m_bfield);
00109
00110
00111
00112 IMdcCalibFunSvc* imdcCalibSvc;
00113 sc = service ("MdcCalibFunSvc", imdcCalibSvc);
00114 m_mdcCalibFunSvc = dynamic_cast<MdcCalibFunSvc*>(imdcCalibSvc);
00115 if ( sc.isFailure() ){
00116 log << MSG::FATAL << "Could not load MdcCalibFunSvc!" << endreq;
00117 return StatusCode::FAILURE;
00118 }
00119
00120
00121 IRawDataProviderSvc* iRawDataProvider;
00122 sc = service ("RawDataProviderSvc", iRawDataProvider);
00123 if ( sc.isFailure() ){
00124 log << MSG::FATAL << "Could not load RawDataProviderSvc!" << endreq;
00125 return StatusCode::FAILURE;
00126 }
00127 m_rawDataProviderSvc = dynamic_cast<RawDataProviderSvc*>(iRawDataProvider);
00128
00129
00130 IMdcUtilitySvc * imdcUtilitySvc;
00131 sc = service ("MdcUtilitySvc", imdcUtilitySvc);
00132 m_mdcUtilitySvc = dynamic_cast<MdcUtilitySvc*> (imdcUtilitySvc);
00133 if ( sc.isFailure() ){
00134 log << MSG::FATAL << "Could not load MdcUtilitySvc!" << endreq;
00135 return StatusCode::FAILURE;
00136 }
00137
00138
00139 IMdcPrintSvc* imdcPrintSvc;
00140 sc = service ("MdcPrintSvc", imdcPrintSvc);
00141 m_mdcPrintSvc = dynamic_cast<MdcPrintSvc*> (imdcPrintSvc);
00142 if ( sc.isFailure() ){
00143 log << MSG::FATAL << "Could not load MdcPrintSvc!" << endreq;
00144 return StatusCode::FAILURE;
00145 }
00146
00147 if(m_hist){
00148
00149 NTuplePtr ntCsmcSew(ntupleSvc(), "MdcxReco/csmc");
00150 if ( ntCsmcSew ) { m_xtupleCsmcSew = ntCsmcSew;}
00151 else {
00152 m_xtupleCsmcSew = ntupleSvc()->book ("MdcxReco/csmc", CLID_ColumnWiseTuple, "MdcxReco reconsturction results");
00153 if ( m_xtupleCsmcSew ) {
00154 sc = m_xtupleCsmcSew->addItem ("dD0", m_csmcD0);
00155 sc = m_xtupleCsmcSew->addItem ("dPhi0", m_csmcPhi0);
00156 sc = m_xtupleCsmcSew->addItem ("dZ0", m_csmcZ0);
00157 sc = m_xtupleCsmcSew->addItem ("dOmega", m_csmcOmega);
00158 sc = m_xtupleCsmcSew->addItem ("dPt", m_csmcPt);
00159 sc = m_xtupleCsmcSew->addItem ("dTanl", m_csmcTanl);
00160 }else{
00161 log << MSG::FATAL << "Could not book MdcxReco/csmc!" << endreq;
00162 return StatusCode::FAILURE;
00163 }
00164 }
00165 }
00166
00167 return StatusCode::SUCCESS;
00168 }
00169
00170 StatusCode MdcxCosmicSewer::execute() {
00171 MsgStream log(msgSvc(), name());
00172 log << MSG::INFO << "in execute()" << endreq;
00173
00174
00175 setFilterPassed(false);
00176
00177 SmartDataPtr<Event::EventHeader> evtHead(eventSvc(),"/Event/EventHeader");
00178 if (!evtHead) {
00179 log << MSG::FATAL<< "Could not retrieve event header" << endreq;
00180 return StatusCode::FAILURE;
00181 }
00182 long t_runNo = evtHead->runNumber();
00183 long t_evtNo = evtHead->eventNumber();
00184 if (m_debug) std::cout << "sew "<<++i_evt<<" run "<<t_runNo<<" evt " << t_evtNo << std::endl;
00185
00186
00187
00188 m_bunchT0 = -999.;
00189 SmartDataPtr<RecEsTimeCol> aevtimeCol(eventSvc(),"/Event/Recon/RecEsTimeCol");
00190 if (!aevtimeCol || aevtimeCol->size()==0) {
00191 log << MSG::WARNING<< "Could not find RecEsTimeCol"<< endreq;
00192 return StatusCode::SUCCESS;
00193 }
00194 RecEsTimeCol::iterator iter_evt = aevtimeCol->begin();
00195 for(; iter_evt!=aevtimeCol->end(); iter_evt++){
00196 m_bunchT0 = (*iter_evt)->getTest();
00197 }
00198
00199
00200 SmartDataPtr<RecMdcTrackCol> recMdcTrackCol(eventSvc(),EventModel::Recon::RecMdcTrackCol);
00201 SmartDataPtr<RecMdcHitCol> recMdcHitCol(eventSvc(),EventModel::Recon::RecMdcHitCol);
00202 if (!recMdcTrackCol || ! recMdcHitCol) return StatusCode::SUCCESS;
00203 if (2!=recMdcTrackCol->size()){
00204 clearRecMdcTrackHit();
00205 return StatusCode::SUCCESS;
00206 }
00207 if(m_debug)std::cout<<name()<<" nTk==2 begin sewing"<< std::endl;
00208
00209
00210
00211 double dParam[5]={0.,0.,0.,0.,0.};
00212 float bprob =0.;
00213 float chisum =0.;
00214 RecMdcTrackCol::iterator besthel;
00215
00216 int besthelId = -999;
00217 RecMdcTrackCol::iterator it = recMdcTrackCol->begin();
00218 for (int iTk=0; it != recMdcTrackCol->end(); it++,iTk++) {
00219 RecMdcTrack* tk = *it;
00220 double Bz = m_bfield->bFieldNominal();
00221 double d0 = -tk->helix(0);
00222 double phi0 = tk->helix(1)+ CLHEP::halfpi;
00223 double omega = Bz*tk->helix(2)/-333.567;
00224 double z0 = tk->helix(3);
00225 double tanl = tk->helix(4);
00226
00227 if(phi0>Constants::pi)phi0-=Constants::twoPi;
00228 if(phi0<-Constants::pi)phi0+=Constants::twoPi;
00229
00230 if(m_debug)std::cout<<__FILE__<<" "<<__LINE__<< "tk"<<iTk
00231 <<"("<<d0<<","<<phi0<<"," <<","<<omega<<","<<z0<<","<<tanl<<")"<<std::endl;
00232
00233 if(iTk==0){
00234 dParam[0]=d0;
00235 dParam[1]=phi0;
00236 dParam[2]=omega;
00237 dParam[3]=z0;
00238 dParam[4]=tanl;
00239 }else{
00240 dParam[0] += d0;
00241 dParam[1] -= (phi0+Constants::pi);
00242 dParam[2] -= omega;
00243 dParam[3] -= z0;
00244 dParam[4] += tanl;
00245 }
00246
00247 if(phi0<0) {
00248 besthelId = iTk;
00249 besthel = it;
00250 }
00251 float bchisq = tk->chi2();
00252 int bndof = tk->ndof();
00253 bprob=Mdcxprobab(bndof,bchisq);
00254 chisum += bchisq;
00255
00256
00257
00258
00259 }
00260
00261 if(besthelId == -999) return StatusCode::SUCCESS;
00262
00263
00264 TrkHelixMaker helixfactory;
00265 TrkLineMaker linefactory;
00266
00267 if(m_debug){
00268 std::cout<<__FILE__<<" param diff: " << "\n D0 " << dParam[0]
00269 << "\n Phi0 " << dParam[1] << "\n Omega " << dParam[2]
00270 << "\n Z0 " << dParam[3] << "\n Tanl " << dParam[4] << std::endl;
00271 }
00272
00273 if(m_hist){
00274 m_csmcD0= dParam[0];
00275 m_csmcPhi0=dParam[1];
00276 m_csmcOmega=dParam[2];
00277 m_csmcZ0=dParam[3];
00278 m_csmcTanl=dParam[4];
00279 m_xtupleCsmcSew->write();
00280 }
00281
00282
00283 if( (fabs(dParam[0])>m_cosmicSewPar[0]) ||
00284 (fabs(dParam[1])>m_cosmicSewPar[1]) ||
00285 (fabs(dParam[2])>m_cosmicSewPar[2]) ||
00286 (fabs(dParam[3])>m_cosmicSewPar[3]) ||
00287 (fabs(dParam[4])>m_cosmicSewPar[4]) ){
00288 if(m_debug)std::cout<<__FILE__<<" 2 trk param not satisfied skip!"<<std::endl;
00289 if(m_debug){
00290 if(fabs(dParam[0])>m_cosmicSewPar[0]) std::cout<<" cut by d0 "<<std::endl;
00291 if(fabs(dParam[1])>m_cosmicSewPar[1]) std::cout<<" cut by phi0 "<<std::endl;
00292 if(fabs(dParam[2])>m_cosmicSewPar[2]) std::cout<<" cut by omega "<<std::endl;
00293 if(fabs(dParam[3])>m_cosmicSewPar[3]) std::cout<<" cut by z0"<<std::endl;
00294 if(fabs(dParam[4])>m_cosmicSewPar[4]) std::cout<<" cut by tanl"<<std::endl;
00295 }
00296 clearRecMdcTrackHit();
00297 return StatusCode::SUCCESS;
00298 }
00299
00300
00301 RecMdcTrack* tk = *besthel;
00302 double Bz = m_bfield->bFieldNominal();
00303 double d0 = -tk->helix(0);
00304 double phi0 = tk->helix(1)+ CLHEP::halfpi;
00305 double omega = Bz*tk->helix(2)/-333.567;
00306 double z0 = tk->helix(3);
00307 double tanl = tk->helix(4);
00308
00309 TrkExchangePar tt(d0,phi0,omega,z0,tanl);
00310 if(m_debug)std::cout<<__FILE__<<" best helix: No."<<tk->trackId()<<" Pat param=("<<d0<<" "<<phi0
00311 <<" "<<omega<<" "<<z0<<" "<<tanl<<")"<< std::endl;
00312
00313
00314
00315 TrkRecoTrk* newTrack;
00316 if(m_lineFit){
00317 newTrack = linefactory.makeTrack(tt, chisum, *m_context, m_bunchT0*1.e-9);
00318 linefactory.setFlipAndDrop(*newTrack, true, true);
00319 }else{
00320 newTrack = helixfactory.makeTrack(tt, chisum, *m_context, m_bunchT0*1.e-9);
00321 helixfactory.setFlipAndDrop(*newTrack, true, true);
00322 }
00323
00324
00325
00326 HitRefVec skippedHits;
00327 it = recMdcTrackCol->begin();
00328 for (int iTk=0; it != recMdcTrackCol->end(); it++,iTk++) {
00329 RecMdcTrack* tk = *it;
00330 HitRefVec hits = tk->getVecHits();
00331 HepVector helixPatPar = m_mdcUtilitySvc->besPar2PatPar(tk->helix());
00332 MdcxHitsToHots(helixPatPar,newTrack,hits,skippedHits);
00333 }
00334
00335
00336
00337
00338
00339 TrkErrCode err=newTrack->hits()->fit();
00340 const TrkFit* newFit = newTrack->fitResult();
00341
00342
00343 if (!newFit) {
00344
00345 if(m_debug)std::cout<<__FILE__<<" sew fit failed!!!"<< std::endl;
00346 HitRefVec::iterator it= skippedHits.begin();
00347 for (;it!=skippedHits.end();++it){ delete it->data();}
00348 delete newTrack;
00349 clearRecMdcTrackHit();
00350 } else {
00351
00352
00353
00354 if(m_lineFit){
00355 linefactory.setFlipAndDrop(*newTrack, false, false);
00356 }else{
00357 helixfactory.setFlipAndDrop(*newTrack, false, false);
00358 }
00359 if(m_debug){
00360 std::cout<<__FILE__<<" sew cosmic fit good "<< std::endl;
00361 newTrack->print(std::cout);
00362 }
00363
00364 SmartDataPtr<MdcHitCol> m_hitCol(eventSvc(), "/Event/Hit/MdcHitCol");
00365 if (!m_hitCol){
00366 m_hitCol= new MdcHitCol;
00367 StatusCode sc = eventSvc()->registerObject("/Event/Hit/MdcHitCol",m_hitCol);
00368 if(!sc.isSuccess()) {
00369 std::cout<< " Could not register MdcHitCol" <<endreq;
00370 return StatusCode::FAILURE;
00371 }
00372 }
00373 uint32_t getDigiFlag = 0;
00374 MdcDigiVec mdcDigiVec = m_rawDataProviderSvc->getMdcDigiVec(getDigiFlag);
00375 const MdcDigi* m_digiMap[43][288];
00376 MdcDigiVec::iterator iter = mdcDigiVec.begin();
00377 for (; iter != mdcDigiVec.end(); iter++ ) {
00378 const MdcDigi* aDigi = *iter;
00379 const Identifier id= aDigi->identify();
00380 int layer = MdcID::layer(id);
00381 int wire = MdcID::wire(id);
00382 m_digiMap[layer][wire] = aDigi;
00383 }
00384
00385
00386 HitRefVec::iterator itHitSkipped = skippedHits.begin();
00387 for (;itHitSkipped!=skippedHits.end();++itHitSkipped){
00388 RecMdcHit* h = *itHitSkipped;
00389 Identifier id(h->getMdcId());
00390 int layer = MdcID::layer(id);
00391 int wire = MdcID::wire(id);
00392 double fltLen = 0.;
00393 HepVector helix = newFit->helix(fltLen).params();
00394 HepSymMatrix err = newFit->helix(fltLen).covariance();
00395 double docaFltLen = h->getFltLen();
00396
00397 double doca = m_mdcUtilitySvc->docaPatPar(layer,wire,helix,err);
00398 double newDoca = fabs(doca);
00399 int flagLR = h->getFlagLR();
00400 if ( flagLR == 0 ){ newDoca = -fabs(doca); }
00401 h->setDoca(newDoca);
00402 if(m_debug>=3)std::cout<<"("<<layer<<","<<wire<<") new doca "<<newDoca<<" resid "<<fabs(fabs(newDoca)-fabs(h->getDriftDistLeft()))<<std::endl;
00403 if(m_debug>=3&&fabs(fabs(newDoca)-fabs(h->getDriftDistLeft()))>0.02) std::cout<<" bad resid "<<fabs(fabs(newDoca)-fabs(h->getDriftDistLeft()))<<" doca "<<newDoca<<" dd "<<h->getDriftDistLeft()<<std::endl;
00404
00405
00406 MdcHit *thehit = new MdcHit(m_digiMap[layer][wire], m_gm);
00407 thehit->setCosmicFit(true);
00408 thehit->setCalibSvc(m_mdcCalibFunSvc);
00409 thehit->setCountPropTime(m_countPropTime);
00410 m_hitCol->push_back(thehit);
00411
00412 HepPoint3D poca; Hep3Vector tempDir;
00413 getInfo(helix,0,poca,tempDir);
00414 if(m_debug>3)std::cout<<"track poca "<<poca<<std::endl;
00415 HepPoint3D pos; Hep3Vector dir;
00416 getInfo(helix,docaFltLen,pos,dir);
00417
00418 if(pos.y()>poca.y()){
00419 int wireAmb = patAmbig(flagLR);
00420
00421 double tof = docaFltLen/Constants::c + (m_bunchT0)*1.e-9;
00422
00423
00424 double eAngle = EntranceAngle(dir.phi() - thehit->phi(pos.z()));
00425 double dAngle = Constants::pi/2 - 1.*dir.theta();
00426 double z = pos.z();
00427 double dist = thehit->driftDist(tof, wireAmb, eAngle, dAngle, z);
00428 double sigma = thehit->sigma(dist,wireAmb,eAngle,dAngle,z);
00429
00430 if(m_debug>3&&fabs(fabs(h->getDoca())-fabs(dist))>0.02)
00431 std::cout<<"("<<layer<<","<<wire<<") old doca "<<h->getDoca()<<" old ddr "<<h->getDriftDistRight()
00432 <<" new dd "<<dist <<" newAmbig "<<wireAmb
00433
00434 <<" old sigma "<<h->getErrDriftDistLeft()
00435 <<" new sigma "<<sigma<<" "<<std::endl;
00436 h->setDriftDistLeft(-1.*abs(dist));
00437 h->setDriftDistRight(abs(dist));
00438 h->setErrDriftDistLeft(sigma);
00439 h->setErrDriftDistRight(sigma);
00440 }
00441 }
00442
00443
00444
00445
00446 store(newTrack, skippedHits);
00447
00448 setFilterPassed(true);
00449 m_nSewed++;
00450 }
00451
00452 if(m_debug) { m_mdcPrintSvc->printRecMdcTrackCol(); }
00453
00454 return StatusCode::SUCCESS;
00455 }
00456
00457 StatusCode MdcxCosmicSewer::finalize() {
00458 MsgStream log(msgSvc(), name());
00459 log << MSG::INFO << "in finalize()" << endreq;
00460 std::cout<<name()<<" after sewed got "<<m_nSewed<<" cosmic tracks"<<std::endl;
00461 return StatusCode::SUCCESS;
00462 }
00463
00464 void MdcxCosmicSewer::MdcxHitsToHots(HepVector& helix, TrkRecoTrk* newTrack, HitRefVec& recMdcHits, HitRefVec& skippedHits) {
00465
00466 TrkHitList* trkHitList = newTrack->hits();
00467
00468 SmartDataPtr<MdcHitCol> m_hitCol(eventSvc(), "/Event/Hit/MdcHitCol");
00469 if (!m_hitCol){
00470 m_hitCol= new MdcHitCol;
00471 StatusCode sc = eventSvc()->registerObject("/Event/Hit/MdcHitCol",m_hitCol);
00472 if(!sc.isSuccess()) {
00473 std::cout<< " Could not register MdcHitCol" <<endreq;
00474 return;
00475 }
00476 }
00477
00478
00479 uint32_t getDigiFlag = 0;
00480 MdcDigiVec mdcDigiVec = m_rawDataProviderSvc->getMdcDigiVec(getDigiFlag);
00481 if (0 == mdcDigiVec.size()){
00482 std::cout << " No hits in MdcDigiVec" << std::endl;
00483 return;
00484 }
00485 const MdcDigi* m_digiMap[43][288];
00486 MdcDigiVec::iterator iter = mdcDigiVec.begin();
00487 for (; iter != mdcDigiVec.end(); iter++ ) {
00488 const MdcDigi* aDigi = *iter;
00489 const Identifier id= aDigi->identify();
00490 int layer = MdcID::layer(id);
00491 int wire = MdcID::wire(id);
00492 m_digiMap[layer][wire] = aDigi;
00493 }
00494
00495
00496 HepPoint3D poca; Hep3Vector tempDir;
00497 getInfo(helix,0,poca,tempDir);
00498 if(m_debug>3)std::cout<<"track poca "<<poca<<std::endl;
00499
00500 HitRefVec::iterator itHit = recMdcHits.begin();
00501 for (;itHit!=recMdcHits.end();++itHit){
00502 RecMdcHit* h = *itHit;
00503 int layer = MdcID::layer(h->getMdcId());
00504 int wire = MdcID::wire(h->getMdcId());
00505
00506 double newFltLen = h->getFltLen();
00507 HepPoint3D pos; Hep3Vector dir;
00508 getInfo(helix,h->getFltLen(),pos,dir);
00509 int newAmbig = patAmbig(h->getFlagLR());
00510 if(pos.y()>poca.y()){
00511 newFltLen *= -1.;
00512 newAmbig *= -1;
00513 if(m_debug>3) std::cout<<" Up track, change sign of fltLen "<<std::endl;
00514 }
00515 int newFlagLR = bes3FlagLR(newAmbig);
00516
00517 if(m_cosmicSewSkip && layer<20){
00518 RecMdcHit* newSkippedHit = new RecMdcHit();
00519 newSkippedHit->setDriftDistLeft(h->getDriftDistLeft());
00520 newSkippedHit->setDriftDistRight(h->getDriftDistRight());
00521 newSkippedHit->setErrDriftDistLeft(h->getErrDriftDistLeft());
00522 newSkippedHit->setErrDriftDistRight(h->getErrDriftDistRight());
00523 newSkippedHit->setChisqAdd(h->getChisqAdd());
00524 newSkippedHit->setFlagLR(newFlagLR);
00525 newSkippedHit->setStat(h->getStat());
00526 newSkippedHit->setMdcId(h->getMdcId());
00527 newSkippedHit->setTdc(h->getTdc());
00528 newSkippedHit->setAdc(h->getAdc());
00529 newSkippedHit->setDriftT(h->getDriftT());
00530 newSkippedHit->setDoca(999.);
00531
00532 newSkippedHit->setEntra(h->getEntra());
00533 newSkippedHit->setZhit(h->getZhit());
00534 newSkippedHit->setFltLen(newFltLen);
00535 skippedHits.push_back(newSkippedHit);
00536 }else{
00537 MdcHit *thehit = new MdcHit(m_digiMap[layer][wire], m_gm);
00538 thehit->setCosmicFit(true);
00539 thehit->setCalibSvc(m_mdcCalibFunSvc);
00540 thehit->setCountPropTime(m_countPropTime);
00541 m_hitCol->push_back(thehit);
00542
00543 MdcRecoHitOnTrack temp(*thehit, newAmbig, m_bunchT0);
00544 MdcHitOnTrack* newhot = &temp;
00545 newhot->setFltLen(newFltLen);
00546 if(m_debug>3) std::cout<<name()<<" ("<<layer<<","<<wire<<") fltLen "<<h->getFltLen()<<" newFltLen "<<newFltLen<<" newAmbig "<<newAmbig<<" pos.y() "<<pos.y()<<std::endl;
00547
00548 trkHitList->appendHot(newhot);
00549 }
00550 }
00551 }
00552
00553 void MdcxCosmicSewer::store(TrkRecoTrk* aTrack, HitRefVec& skippedHits) {
00554 StatusCode sc;
00555 MsgStream log(msgSvc(), name());
00556
00557 assert (aTrack != NULL);
00558
00559
00560 SmartIF<IDataManagerSvc> dataManSvc(eventSvc());
00561 DataObject *aTrackCol;
00562 eventSvc()->findObject("/Event/Recon/RecMdcTrackCol",aTrackCol);
00563 if(aTrackCol != NULL) {
00564 dataManSvc->clearSubTree("/Event/Recon/RecMdcTrackCol");
00565 eventSvc()->unregisterObject("/Event/Recon/RecMdcTrackCol");
00566 }
00567
00568 DataObject *aRecHitCol;
00569 eventSvc()->findObject("/Event/Recon/RecMdcHitCol",aRecHitCol);
00570 if(aRecHitCol != NULL) {
00571 dataManSvc->clearSubTree("/Event/Recon/RecMdcHitCol");
00572 eventSvc()->unregisterObject("/Event/Recon/RecMdcHitCol");
00573 }
00574
00575 SmartDataPtr<RecMdcTrackCol> trackList(eventSvc(),EventModel::Recon::RecMdcTrackCol);
00576 if (!trackList) {
00577 RecMdcTrackCol *trackList = new RecMdcTrackCol;
00578 sc = eventSvc()->registerObject(EventModel::Recon::RecMdcTrackCol, trackList);
00579 if(!sc.isSuccess()) {
00580 log << MSG::FATAL << " Could not register RecMdcTrack collection" <<endreq;
00581 }
00582 }
00583
00584 SmartDataPtr<RecMdcHitCol> hitList(eventSvc(),EventModel::Recon::RecMdcHitCol);
00585 if (!hitList) {
00586 hitList = new RecMdcHitCol;
00587 sc = eventSvc()->registerObject(EventModel::Recon::RecMdcHitCol, hitList);
00588 if(!sc.isSuccess()) {
00589 log << MSG::FATAL << " Could not register RecMdcHit collection" <<endreq;
00590 }
00591 }
00592 int trackId = trackList->size() - 1;
00593
00594 if(m_debug)std::cout<<__FILE__<<" "<<__LINE__<<" STORED"<< std::endl;
00595 MdcTrack mdcTrack(aTrack);
00596
00597
00598 int tkStat = -1;
00599 mdcTrack.storeTrack(trackId, trackList, hitList, tkStat);
00600
00601
00602 RecMdcTrackCol::iterator it = trackList->begin();
00603 HitRefVec hl = (*it)->getVecHits();
00604 HitRefVec hitRefVec;
00605 HitRefVec::iterator itHit = hl.begin();
00606 int ihl=0;
00607 for (;itHit!=hl.end();++itHit){
00608 ihl++;
00609 RecMdcHit* h = *itHit;
00610 SmartRef<RecMdcHit> refHit(h);
00611 hitRefVec.push_back(refHit);
00612 }
00613
00614 if(m_debug) std::cout<<__FILE__<<" "<<__LINE__<<"refHit stored "<< ihl<< std::endl;
00615
00616 HitRefVec::iterator itHitSkipped = skippedHits.begin();
00617 int iskipped=0;
00618 for (;itHitSkipped!=skippedHits.end();++itHitSkipped){
00619 iskipped++;
00620 (*itHitSkipped)->setTrkId(trackId);
00621 hitList->push_back(*itHitSkipped);
00622 SmartRef<RecMdcHit> refHitSkipped(*itHitSkipped);
00623 hitRefVec.push_back(refHitSkipped);
00624 }
00625 if(m_debug) std::cout<<__FILE__<<" "<<__LINE__<<"skippedHits stored "<< iskipped<< std::endl;
00626 (*it)->setVecHits(hitRefVec);
00627
00628 (*it)->setNhits((*it)->getVecHits().size());
00629 (*it)->setNster(-1);
00630
00631
00632 if((*it)->helix(2)<0.00001)(*it)->setNdof((*it)->getNhits()-4);
00633 else (*it)->setNdof((*it)->getNhits()-5);
00634
00635 if(m_debug) std::cout<<__FILE__<<" "<<__LINE__<<" stored "<< (*it)->getVecHits().size()<< std::endl;
00636
00637 }
00638
00639 void MdcxCosmicSewer::clearRecMdcTrackHit(){
00640 MsgStream log(msgSvc(), name());
00641 StatusCode sc;
00642
00643
00644
00645
00646 SmartIF<IDataManagerSvc> dataManSvc(eventSvc());
00647 DataObject *aTrackCol;
00648 eventSvc()->findObject("/Event/Recon/RecMdcTrackCol",aTrackCol);
00649 if(aTrackCol != NULL) {
00650 dataManSvc->clearSubTree("/Event/Recon/RecMdcTrackCol");
00651 eventSvc()->unregisterObject("/Event/Recon/RecMdcTrackCol");
00652 }
00653 DataObject *aRecHitCol;
00654 eventSvc()->findObject("/Event/Recon/RecMdcHitCol",aRecHitCol);
00655 if(aRecHitCol != NULL) {
00656 dataManSvc->clearSubTree("/Event/Recon/RecMdcHitCol");
00657 eventSvc()->unregisterObject("/Event/Recon/RecMdcHitCol");
00658 }
00659 SmartDataPtr<RecMdcTrackCol> trackListTemp(eventSvc(),EventModel::Recon::RecMdcTrackCol);
00660 if (!trackListTemp) {
00661 RecMdcTrackCol *trackListTemp= new RecMdcTrackCol;
00662 sc = eventSvc()->registerObject(EventModel::Recon::RecMdcTrackCol, trackListTemp);
00663 if(!sc.isSuccess()) {
00664 log << MSG::FATAL << " Could not register RecMdcTrack collection" <<endreq;
00665 }
00666 }
00667 SmartDataPtr<RecMdcHitCol> hitListTemp(eventSvc(),EventModel::Recon::RecMdcHitCol);
00668 if (!hitListTemp) {
00669 RecMdcHitCol *hitListTemp= new RecMdcHitCol;
00670 sc = eventSvc()->registerObject(EventModel::Recon::RecMdcHitCol, hitListTemp);
00671 if(!sc.isSuccess()) {
00672 log << MSG::FATAL << " Could not register RecMdcHit collection" <<endreq;
00673 }
00674 }
00675 }
00676
00677 void MdcxCosmicSewer::getInfo(HepVector helix, double fltLen, HepPoint3D& pos, Hep3Vector & dir){
00678 if(m_lineFit) return;
00679 double d0 = helix[0];
00680 double phi0 = helix[1];
00681 double omega = helix[2];
00682 double z0 = helix[3];
00683 double tanDip = helix[4];
00684 double cDip = 1./sqrt(1.+tanDip*tanDip);
00685 double sDip = tanDip * cDip;
00686 double phi00 = phi0;
00687 double ang = phi00 + cDip*fltLen*omega;
00688 double cang = cos(ang);
00689 double sang = sin(ang);
00690 double sphi0 = sin(phi00);
00691 double cphi0 = cos(phi00);
00692 HepPoint3D referencePoint(0.,0.,0.);
00693 double xt = (sang - sphi0)/omega - d0*sphi0 + referencePoint.x();
00694 double yt = -(cang - cphi0)/omega + d0*cphi0 + referencePoint.y();
00695 double zt = z0 + fltLen*sDip + referencePoint.z();
00696 pos.setX(xt);
00697 pos.setY(yt);
00698 pos.setZ(zt);
00699
00700 dir.setX(cang * cDip);
00701 dir.setY(sang * cDip);
00702 dir.setZ(sDip);
00703 }
00704
00705 int MdcxCosmicSewer::patAmbig(int bes3FlagLR){
00706 int ambig = 0;
00707 if ( bes3FlagLR == 0 ){
00708 ambig = 1;
00709 }else if( bes3FlagLR == 1){
00710 ambig = -1;
00711 }else if( bes3FlagLR == 2){
00712 ambig = 0;
00713 }
00714 return ambig;
00715 }
00716
00717 int MdcxCosmicSewer::bes3FlagLR(int patAmbig){
00718 int flagLR = 2;
00719 if ( patAmbig == 1 ){
00720 flagLR = 0;
00721 }else if( patAmbig == -1){
00722 flagLR = 1;
00723 }else if( patAmbig == 0){
00724 flagLR = 2;
00725 }
00726 return flagLR;
00727 }