00001 #include "GaudiKernel/MsgStream.h"
00002 #include "GaudiKernel/AlgFactory.h"
00003 #include "GaudiKernel/ISvcLocator.h"
00004 #include "GaudiKernel/SmartDataPtr.h"
00005 #include "GaudiKernel/IDataProviderSvc.h"
00006 #include "GaudiKernel/IDataManagerSvc.h"
00007 #include "GaudiKernel/PropertyMgr.h"
00008 #include "GaudiKernel/INTupleSvc.h"
00009 #include "EventModel/EventHeader.h"
00010 #include "ReconEvent/ReconEvent.h"
00011 #include "ExtEvent/RecExtTrack.h"
00012 #include "MdcRecEvent/RecMdcTrack.h"
00013 #include "EmcRecEventModel/RecEmcEventModel.h"
00014 #include "MdcRawEvent/MdcDigi.h"
00015 #include "TofRawEvent/TofDigi.h"
00016 #include "EmcRawEvent/EmcDigi.h"
00017 #include "MucRawEvent/MucDigi.h"
00018 #include "McTruth/McKine.h"
00019 #include "McTruth/McParticle.h"
00020 #include "McTruth/MdcMcHit.h"
00021 #include "McTruth/TofMcHit.h"
00022 #include "McTruth/EmcMcHit.h"
00023 #include "McTruth/MucMcHit.h"
00024 #include "McTruth/McRelTableDefs.h"
00025 #include "Identifier/Identifier.h"
00026 #include "Identifier/MucID.h"
00027 #include "GaudiKernel/IPartPropSvc.h"
00028
00029 #include "MucRecAlg/MucRecTrkExt.h"
00030 #include "MucRecAlg/MucRecTrkExtParameter.h"
00031 #include "MucGeomSvc/IMucGeomSvc.h"
00032 #include "MucGeomSvc/MucGeoGeneral.h"
00033 #include "MucGeomSvc/MucGeoGap.h"
00034 #include "MucGeomSvc/MucGeoStrip.h"
00035 #include "MucGeomSvc/MucConstant.h"
00036 #include "MucRecEvent/MucRecHit.h"
00037 #include "MucRecEvent/MucRecHitContainer.h"
00038 #include "MucRecEvent/RecMucTrack.h"
00039 #include <vector>
00040 #include <iostream>
00041 #include <fstream>
00042 #include <iomanip>
00043
00044 using namespace std;
00045 using namespace Event;
00046
00048
00049 MucRecTrkExt::MucRecTrkExt(const std::string& name, ISvcLocator* pSvcLocator) :
00050 Algorithm(name, pSvcLocator)
00051 {
00052
00053 declareProperty("ExtTrackSeedMode", m_ExtTrackSeedMode = 1);
00054 declareProperty("CompareWithMcHit", m_CompareWithMcHit = 0);
00055 declareProperty("FittingMethod", m_fittingMethod = 1);
00056 declareProperty("EmcShowerSeedMode",m_EmcShowerSeedMode = 0);
00057 declareProperty("MucHitSeedMode", m_MucHitSeedMode = 0);
00058 declareProperty("ConfigFile", m_configFile = "MucConfig.xml");
00059 declareProperty("Blind", m_Blind = false);
00060 declareProperty("NtOutput", m_NtOutput = 0);
00061 declareProperty("MsOutput", m_MsOutput = false);
00062 declareProperty("FilterFile", m_filter_filename);
00063
00064 }
00065
00066
00067 StatusCode MucRecTrkExt::initialize()
00068 {
00069 MsgStream log(msgSvc(), name());
00070 log << MSG::INFO << "in initialize()" << endreq;
00071
00072
00073 IPartPropSvc* p_PartPropSvc;
00074 static const bool CREATEIFNOTTHERE(true);
00075 StatusCode PartPropStatus = service("PartPropSvc", p_PartPropSvc, CREATEIFNOTTHERE);
00076 if (!PartPropStatus.isSuccess() || 0 == p_PartPropSvc) {
00077 log << MSG::WARNING << " Could not initialize Particle Properties Service" << endreq;
00078 return PartPropStatus;
00079 }
00080 m_particleTable = p_PartPropSvc->PDT();
00081
00082 m_totalEvent = 0;
00083
00084 m_NDigisTotal = 0;
00085 m_NHitsTotal = 0;
00086 m_NHitsFoundTotal = 0;
00087 m_NHitsLostTotal = 0;
00088 m_NHitsMisFoundTotal = 0;
00089 m_NHitsLostByMdcTotal = 0;
00090 m_NHitsLostByExtTotal = 0;
00091
00092 m_NTracksTotal = 0;
00093 m_NTracksRecoTotal = 0;
00094 m_NTracksLostByMdcTotal = 0;
00095 m_NTracksLostByExtTotal = 0;
00096 m_NTracksMdcGhostTotal = 0;
00097
00098 for(int i = 0; i < 20; i++) m_NHitsLost.push_back(0);
00099 for(int i = 0; i < 10; i++) m_NHitsLostInGap.push_back(0);
00100
00101 IMucGeomSvc* mucGeomSvc;
00102 StatusCode sc = service("MucGeomSvc", mucGeomSvc);
00103 if (sc == StatusCode::SUCCESS) {
00104
00105 mucGeomSvc->Dump();
00106
00107 } else {
00108 return StatusCode::FAILURE;
00109 }
00110 m_MucRecHitContainer = 0;
00111
00112
00113
00114 if(m_NtOutput>=1){
00115 NTuplePtr nt1(ntupleSvc(), "FILE401/T");
00116
00117 if ( nt1 ) { m_tuple = nt1;}
00118 else {
00119
00120 m_tuple = ntupleSvc()->book ("FILE401/T", CLID_RowWiseTuple, "MucTrkRecon N-Tuple");
00121 if ( m_tuple ) {
00122 sc = m_tuple->addItem ("posx", m_posx);
00123 sc = m_tuple->addItem ("posy", m_posy);
00124 sc = m_tuple->addItem ("posz", m_posz);
00125 sc = m_tuple->addItem ("posx_ext", m_posx_ext);
00126 sc = m_tuple->addItem ("posy_ext", m_posy_ext);
00127 sc = m_tuple->addItem ("posz_ext", m_posz_ext);
00128 sc = m_tuple->addItem ("posxsigma", m_posx_sigma);
00129 sc = m_tuple->addItem ("posysigma", m_posy_sigma);
00130 sc = m_tuple->addItem ("poszsigma", m_posz_sigma);
00131 sc = m_tuple->addItem ("Depth", m_depth);
00132 sc = m_tuple->addItem ("Distance",m_Distance);
00133 sc = m_tuple->addItem ("MaxHits", m_MaxHits);
00134 sc = m_tuple->addItem ("Chi2", m_Chi2);
00135 sc = m_tuple->addItem ("Dist_x", m_Dist_x);
00136 sc = m_tuple->addItem ("Dist_y", m_Dist_y);
00137 sc = m_tuple->addItem ("Dist_z", m_Dist_z);
00138 sc = m_tuple->addItem ("px_mc", m_px_mc);
00139 sc = m_tuple->addItem ("py_mc", m_py_mc);
00140 sc = m_tuple->addItem ("pz_mc", m_pz_mc);
00141 sc = m_tuple->addItem ("emctrack",m_emctrack);
00142 sc = m_tuple->addItem ("muchasdigi",m_muc_digi);
00143
00144 sc = m_tuple->addItem ("part", m_part);
00145 sc = m_tuple->addItem ("seg", m_seg);
00146 sc = m_tuple->addItem ("gap", m_gap);
00147 sc = m_tuple->addItem ("strip", m_strip);
00148 sc = m_tuple->addItem ("diff", m_diff);
00149 sc = m_tuple->addItem ("distance",m_distance);
00150 sc = m_tuple->addItem ("run", m_run);
00151 sc = m_tuple->addItem ("event", m_event);
00152
00153 }
00154 else {
00155 log << MSG::ERROR << " Cannot book N-tuple:" << long(m_tuple) << endmsg;
00156
00157 }
00158 }
00159 }
00160
00161
00162 if ( m_filter_filename == "") {
00163 m_filter_filename =getenv("MUCRECALGROOT");
00164 m_filter_filename +="/share/filter.txt";
00165 }
00166
00167 if (m_filter_filename.size()) {
00168 std::ifstream infile(m_filter_filename.c_str());
00169
00170 while (!infile.eof()) {
00171 FilterEvent filterevt;
00172 infile >> filterevt.bossver >> filterevt.runid >> filterevt.eventid;
00173 if ((!infile.good()) || (infile.eof())) {
00174 break;
00175 }
00176
00177 m_filter_event.push_back(filterevt);
00178
00179
00180
00181 }
00182 }
00183
00184
00185
00186 return StatusCode::SUCCESS;
00187 }
00188
00189
00190 StatusCode MucRecTrkExt::execute() {
00191
00192 MsgStream log(msgSvc(), name());
00193 log << MSG::INFO << "in execute()" << endreq;
00194
00195
00196 StatusCode sc = StatusCode::SUCCESS;
00197
00198
00199 int event, run;
00200
00201 SmartDataPtr<Event::EventHeader> eventHeader(eventSvc(),"/Event/EventHeader");
00202 if (!eventHeader) {
00203 log << MSG::FATAL << "Could not find Event Header" << endreq;
00204
00205 }
00206 m_totalEvent ++;
00207 log << MSG::INFO << "Event: "<<m_totalEvent<<"\tcurrent run: "<<eventHeader->runNumber()<<"\tcurrent event: "<< eventHeader->eventNumber()<<endreq;
00208
00209 event = eventHeader->eventNumber();
00210 run = eventHeader->runNumber();
00211
00212 string release = getenv("BES_RELEASE");
00213
00214
00215
00216 for (std::vector<FilterEvent>::iterator it = m_filter_event.begin();
00217 it != m_filter_event.end(); ++it) {
00218 const FilterEvent& fe = (*it);
00219 if (release == fe.bossver && run == fe.runid && event == fe.eventid) {
00220 cout << "SKIP: " << fe.bossver << " "
00221 << fe.runid << " "
00222 << fe.eventid << std::endl;
00223 return StatusCode::SUCCESS;
00224 }
00225 }
00226
00227
00228 if(m_CompareWithMcHit==1){
00229 SmartDataPtr<McParticleCol> mcParticleCol(eventSvc(),"/Event/MC/McParticleCol");
00230 if (!mcParticleCol) {
00231 log << MSG::FATAL << "Could not find McParticle" << endreq;
00232
00233 }
00234
00235 McParticleCol::iterator iter_mc = mcParticleCol->begin();
00236
00237
00238
00239 int pid = (*iter_mc)->particleProperty();
00240 int charge = 0;
00241 double mass = -1.0;
00242
00243 if( pid >0 )
00244 {
00245 if(m_particleTable->particle( pid )){
00246 charge = (int)m_particleTable->particle( pid )->charge();
00247 mass = m_particleTable->particle( pid )->mass();
00248 }
00249 else log << MSG::ERROR << "no this particle id in particle table, please check data" <<endreq;
00250 } else if ( pid <0 ) {
00251 if(m_particleTable->particle( -pid )){
00252 charge = (int)m_particleTable->particle( -pid )->charge();
00253 charge *= -1;
00254 mass = m_particleTable->particle( -pid )->mass();
00255 }
00256 else log << MSG::ERROR << "no this particle id in particle table, please check data" <<endreq;
00257 } else {
00258 log << MSG::WARNING << "wrong particle id, please check data" <<endreq;
00259 }
00260
00261
00262
00263
00264
00265
00266
00267 HepLorentzVector initialMomentum = (*iter_mc)->initialFourMomentum();
00268 HepLorentzVector initialPos = (*iter_mc)->initialPosition();
00269 if(m_NtOutput>=1){
00270 m_px_mc = initialMomentum.px();
00271 m_py_mc = initialMomentum.py();
00272 m_pz_mc = initialMomentum.pz();
00273 }
00274
00275
00276 log << MSG::INFO << " particleId = " << (*iter_mc)->particleProperty() << endreq;
00277 }
00278
00279
00280 int digiId = 0;
00281 SmartDataPtr<MucDigiCol> mucDigiCol(eventSvc(),"/Event/Digi/MucDigiCol");
00282 if (!mucDigiCol) {
00283 log << MSG::FATAL << "Could not find MUC digi" << endreq;
00284 return( StatusCode::FAILURE);
00285 }
00286 MucDigiCol::iterator iter4 = mucDigiCol->begin();
00287 for (;iter4 != mucDigiCol->end(); iter4++, digiId++) {
00288 }
00289 log << MSG::INFO << "Total MUC digis:\t" << digiId << endreq;
00290 if( digiId == 0 ) return ( StatusCode::SUCCESS);
00291
00292
00293
00294
00295
00296
00297
00298
00299
00300
00301
00302
00303
00304
00305
00306
00307
00308
00309
00310
00311
00312
00313
00314
00315
00316
00317
00318 Identifier mucID;
00319
00320
00321
00322
00323 if (m_CompareWithMcHit) {
00324 McPartToMucHitTab assocMcMuc;
00325 assocMcMuc.init();
00326
00327 SmartDataPtr<McParticleCol> mcParticleCol(eventSvc(),"/Event/MC/McParticleCol");
00328 if (!mcParticleCol) {
00329 log << MSG::FATAL << "Could not find McParticle" << endreq;
00330
00331 }
00332 McParticleCol::iterator iter_mc = mcParticleCol->begin();
00333
00334
00335 SmartDataPtr<Event::MucMcHitCol> aMucMcHitCol(eventSvc(),"/Event/MC/MucMcHitCol");
00336 if (!aMucMcHitCol) {
00337 log << MSG::WARNING << "Could not find MucMcHitCol" << endreq;
00338
00339 }
00340
00341 log << MSG::DEBUG << "MucMcHitCol contains " << aMucMcHitCol->size() << " Hits " << endreq;
00342 MucMcHitCol::iterator iter_MucMcHit = aMucMcHitCol->begin();
00343 for (; iter_MucMcHit != aMucMcHitCol->end(); iter_MucMcHit++) {
00344 mucID = (*iter_MucMcHit)->identify();
00345
00346 log << MSG::DEBUG << " MucMcHit " << " : "
00347 << " part " << MucID::barrel_ec(mucID)
00348 << " seg " << MucID::segment(mucID)
00349 << " gap " << MucID::layer(mucID)
00350 << " strip " << MucID::channel(mucID)
00351 << " Track Id " << (*iter_MucMcHit)->getTrackIndex()
00352 << " pos x " << (*iter_MucMcHit)->getPositionX()
00353 << " pos y " << (*iter_MucMcHit)->getPositionY()
00354 << " pos z " << (*iter_MucMcHit)->getPositionZ()
00355 << endreq;
00356
00357 McParticle *assocMcPart = 0;
00358 for (iter_mc = mcParticleCol->begin(); iter_mc != mcParticleCol->end(); iter_mc++) {
00359 if ( (*iter_mc)->trackIndex() == (*iter_MucMcHit)->getTrackIndex() ) {
00360 assocMcPart = *iter_mc;
00361 break;
00362 }
00363 }
00364 if (assocMcPart == 0) {
00365 log << MSG::WARNING << "No Corresponding Mc Particle found for this MucMcHit" << endreq;
00366 }
00367
00368 MucMcHit *assocMucMcHit = *iter_MucMcHit;
00369 McPartToMucHitRel *relMcMuc = 0;
00370 relMcMuc = new McPartToMucHitRel( assocMcPart, assocMucMcHit );
00371 if (relMcMuc == 0) log << MSG::DEBUG << "relMcMuc not created " << endreq;
00372 else {
00373 bool addstat = assocMcMuc.addRelation( relMcMuc );
00374 if(!addstat) delete relMcMuc;
00375 }
00376 }
00377
00378 log << MSG::DEBUG << " Fill McPartToMucHitTab, size " << assocMcMuc.size() << endreq;
00379 iter_mc = mcParticleCol->begin();
00380 for (;iter_mc != mcParticleCol->end(); iter_mc++) {
00381 log << MSG::DEBUG << " Track index " << (*iter_mc)->trackIndex()
00382 << " particleId = " << (*iter_mc)->particleProperty()
00383 << endreq;
00384
00385 vector< McPartToMucHitRel* > vecMucMcHit = assocMcMuc.getRelByFirst(*iter_mc);
00386 vector< McPartToMucHitRel* >::iterator iter_MucMcHit = vecMucMcHit.begin();
00387 for (; iter_MucMcHit != vecMucMcHit.end(); iter_MucMcHit++) {
00388 mucID = (*iter_MucMcHit)->getSecond()->identify();
00389
00390 log << MSG::DEBUG
00391
00392 << " MC Particle index "
00393 << (*iter_MucMcHit)->getFirst()->trackIndex()
00394 << " contains " << " MucMcHit "
00395 << " part " << MucID::barrel_ec(mucID)
00396 << " seg " << MucID::segment(mucID)
00397 << " gap " << MucID::layer(mucID)
00398 << " strip " << MucID::channel(mucID)
00399
00400
00401
00402 << endreq;
00403 }
00404 }
00405
00406
00407
00408
00409
00410
00411
00412
00413
00414
00415
00416
00417 }
00418
00419
00420
00421 if (!m_MucRecHitContainer) {
00422 m_MucRecHitContainer = new MucRecHitContainer();
00423 }
00424 m_MucRecHitContainer->Clear();
00425 MucRecHitCol *aMucRecHitCol = new MucRecHitCol();
00426 m_MucRecHitContainer->SetMucRecHitCol(aMucRecHitCol);
00427
00428 SmartIF<IDataManagerSvc> dataManSvc(eventSvc());
00429 DataObject *mucRecHitCol;
00430 eventSvc()->findObject("/Event/Recon/MucRecHitCol",mucRecHitCol);
00431 if(mucRecHitCol != NULL) {
00432 dataManSvc->clearSubTree("/Event/Recon/MucRecHitCol");
00433 eventSvc()->unregisterObject("/Event/Recon/MucRecHitCol");
00434 }
00435
00436 sc = eventSvc()->registerObject("/Event/Recon/MucRecHitCol", aMucRecHitCol);
00437 if (!sc) {
00438 log << MSG::ERROR << "/Event/Recon/MucRecHitCol not registerd!" << endreq;
00439 return( StatusCode::FAILURE);
00440 }
00441
00442 log << MSG::INFO << "Add digis" << endreq;
00443 MucDigiCol::iterator iter_Muc = mucDigiCol->begin();
00444 int mucDigiId = 0;
00445 for (;iter_Muc != mucDigiCol->end(); iter_Muc++, mucDigiId++) {
00446 mucID = (*iter_Muc)->identify();
00447 int part = MucID::barrel_ec(mucID);
00448 int seg = MucID::segment(mucID);
00449 int gap = MucID::layer(mucID);
00450 int strip = MucID::channel(mucID);
00451
00452 m_MucRecHitContainer->AddHit(part, seg, gap, strip);
00453
00454 log << MSG::DEBUG << " digit" << mucDigiId << " : "
00455 << " part " << part
00456 << " seg " << seg
00457 << " gap " << gap
00458 << " strip " << strip
00459 << endreq;
00460 }
00461
00462
00463
00464 RecMucTrackCol *aRecMucTrackCol = new RecMucTrackCol();
00465
00466
00467 DataObject *aReconEvent ;
00468 eventSvc()->findObject("/Event/Recon",aReconEvent);
00469 if(aReconEvent==NULL) {
00470
00471 aReconEvent = new ReconEvent();
00472 StatusCode sc = eventSvc()->registerObject("/Event/Recon",aReconEvent);
00473 if(sc!=StatusCode::SUCCESS) {
00474 log << MSG::FATAL << "Could not register ReconEvent" <<endreq;
00475 return( StatusCode::FAILURE);
00476 }
00477 }
00478 StatusCode fr = eventSvc()->findObject("/Event/Recon",aReconEvent);
00479 if(fr!=StatusCode::SUCCESS) {
00480 log << MSG::WARNING << "Could not find register ReconEvent, will create it" <<endreq;
00481 StatusCode sc = eventSvc()->registerObject("/Event/Recon",aReconEvent);
00482 if(sc!=StatusCode::SUCCESS) {
00483 log << MSG::FATAL << "Could not register ReconEvent" <<endreq;
00484 return( StatusCode::FAILURE);
00485 }
00486 }
00487
00488 DataObject *mucTrackCol;
00489 eventSvc()->findObject("/Event/Recon/RecMucTrackCol",mucTrackCol);
00490 if(mucTrackCol != NULL) {
00491 dataManSvc->clearSubTree("/Event/Recon/RecMucTrackCol");
00492 eventSvc()->unregisterObject("/Event/Recon/RecMucTrackCol");
00493 }
00494
00495 sc = eventSvc()->registerObject("/Event/Recon/RecMucTrackCol", aRecMucTrackCol);
00496 if(sc!=StatusCode::SUCCESS) {
00497 log << MSG::FATAL << "Could not register MUC track collection" << endreq;
00498 return( StatusCode::FAILURE);
00499 }
00500
00501
00502 SmartDataPtr<RecMucTrackCol> findRecMucTrackCol(eventSvc(),"/Event/Recon/RecMucTrackCol");
00503 if (!findRecMucTrackCol) {
00504 log << MSG::FATAL << "Could not find RecMucTrackCol" << endreq;
00505 return( StatusCode::FAILURE);
00506 }
00507 aRecMucTrackCol->clear();
00508
00509
00510
00511
00512
00513 SmartDataPtr<RecExtTrackCol> aExtTrackCol(eventSvc(),"/Event/Recon/RecExtTrackCol");
00514 if (!aExtTrackCol) {
00515 log << MSG::WARNING << "Can't find ExtTrackCol in TDS!" << endreq;
00516
00517 }
00518
00519 SmartDataPtr<RecMdcTrackCol> aMdcTrackCol(eventSvc(),"/Event/Recon/RecMdcTrackCol");
00520 if (!aMdcTrackCol) {
00521 log << MSG::WARNING << "Can't find MdcTrackCol in TDS!" << endreq;
00522
00523 }
00524
00525
00526
00527 SmartDataPtr<RecEmcShowerCol> aShowerCol(eventSvc(),"/Event/Recon/RecEmcShowerCol");
00528 if (!aShowerCol) {
00529 log << MSG::WARNING << "Could not find RecEmcShowerCol" << endreq;
00530
00531 }
00532
00533
00534
00535
00536
00537
00538
00539
00540
00541 int muctrackId = 0;
00542
00543 log << MSG::INFO << "Add tracks info, ExtTrackSeedMode = " << m_ExtTrackSeedMode << ", EmcShowerSeedMode = " << m_EmcShowerSeedMode << endreq;
00544
00545 if (m_ExtTrackSeedMode == 1)
00546 {
00547 SmartDataPtr<McParticleCol> mcParticleCol(eventSvc(),"/Event/MC/McParticleCol");
00548 if (!mcParticleCol) {
00549 log << MSG::FATAL << "Could not find McParticle" << endreq;
00550
00551 return( StatusCode::SUCCESS);
00552 }
00553 McParticleCol::iterator iter_mc = mcParticleCol->begin();
00554
00555 int trackIndex = -99;
00556 for (int iTrack = 0;iter_mc != mcParticleCol->end(); iter_mc++, iTrack++) {
00557 if(!(*iter_mc)->primaryParticle()) continue;
00558 int pid = (*iter_mc)->particleProperty();
00559 int charge = 0;
00560 double mass = -1.0;
00561 if( pid >0 )
00562 {
00563 if(m_particleTable->particle( pid )){
00564 charge = (int)m_particleTable->particle( pid )->charge();
00565 mass = m_particleTable->particle( pid )->mass();
00566 }
00567 else log << MSG::ERROR << "no this particle id in particle table, please check data" <<endreq;
00568 } else if ( pid <0 )
00569 {
00570 if(m_particleTable->particle( -pid )){
00571 charge = (int)m_particleTable->particle( -pid )->charge();
00572 charge *= -1;
00573 mass = m_particleTable->particle( -pid )->mass();
00574 }
00575 else log << MSG::ERROR << "no this particle id in particle table, please check data" <<endreq;
00576 } else {
00577 log << MSG::WARNING << "wrong particle id, please check data" <<endreq;
00578 }
00579
00580 if(!charge) {
00581 log << MSG::WARNING
00582 << " neutral particle charge = 0!!! ...just skip it !"
00583 << endreq;
00584 continue;
00585 }
00586
00587 trackIndex = (*iter_mc)->trackIndex();
00588 log << MSG::DEBUG << "iTrack " << iTrack << " index " << trackIndex
00589 << " particleId = " << (*iter_mc)->particleProperty()
00590 << endreq;
00591
00592 RecMucTrack *aTrack = new RecMucTrack();
00593 aTrack->setTrackId(trackIndex);
00594 aTrack->setId(muctrackId);
00595
00596 HepLorentzVector initialMomentum = (*iter_mc)->initialFourMomentum();
00597 HepLorentzVector initialPos = (*iter_mc)->initialPosition();
00598 float theta0 = initialMomentum.theta();
00599 float phi0 = initialMomentum.phi();
00600
00601
00602
00603 float x0 = initialPos.x();
00604 float y0 = initialPos.y();
00605 float z0 = initialPos.z();
00606
00607 Hep3Vector startPos(x0, y0, z0);
00608 Hep3Vector startP(initialMomentum.px(), initialMomentum.py(), initialMomentum.pz());
00609 log << MSG::DEBUG << "startP " << startP <<" startPos "<<startPos<< endreq;
00610 Hep3Vector endPos(0, 0, 0), endP;
00611
00612 aTrack->GetMdcExtTrack(startPos, startP, charge, endPos, endP);
00613 aTrack->SetMdcPos( x0, y0, z0);
00614 aTrack->SetMdcMomentum( startP.x(), startP.y(), startP.z() );
00615 aTrack->SetExtTrackID(trackIndex);
00616 aTrack->SetExtMucPos( endPos.x(), endPos.y(), endPos.z() );
00617 aTrack->SetExtMucMomentum( endP.x(), endP.y(), endP.z() );
00618 aTrack->SetMucPos( endPos.x(), endPos.y(), endPos.z() );
00619 aTrack->SetMucMomentum( endP.x(), endP.y(), endP.z() );
00620 aTrack->SetCurrentPos( endPos.x(), endPos.y(), endPos.z());
00621 aTrack->SetCurrentDir( endP.x(), endP.y(), endP.z() );
00622
00623
00624
00625
00626
00627
00628
00629
00630
00631
00632 log << MSG::INFO << "Add track, trkIndex: " << trackIndex << "\tmucTrkId: " << muctrackId << endreq;
00633 aRecMucTrackCol->add(aTrack);
00634 muctrackId ++ ;
00635 }
00636 }
00637 else if (m_ExtTrackSeedMode == 2)
00638 {
00639 if (!aExtTrackCol || !aMdcTrackCol) {
00640 log << MSG::WARNING << "Can't find ExtTrackCol or MdcTrackCol in TDS!" << endreq;
00641 return StatusCode::SUCCESS;
00642 }
00643
00644 int trackIndex = -99;
00645 double kdep = -99;
00646 double krechi = 0.;
00647 int kdof = 0;
00648 int kbrLay = -1;
00649 int kecLay = -1;
00650
00651 RecExtTrackCol::iterator iter_ExtTrack = aExtTrackCol->begin();
00652 RecMdcTrackCol::iterator iter_MdcTrack = aMdcTrackCol->begin();
00653 int iExtTrack = 0;
00654 for (;iter_ExtTrack != aExtTrackCol->end() && iter_MdcTrack != aMdcTrackCol->end(); iter_ExtTrack++, iter_MdcTrack++, iExtTrack++)
00655 {
00656 trackIndex = (*iter_ExtTrack)->GetTrackId();
00657 log << MSG::DEBUG << "iExtTrack " << iExtTrack << " index " << trackIndex << " MucPos "
00658 << iExtTrack << (*iter_ExtTrack)->mucPosition().x() << " "
00659 << (*iter_ExtTrack)->mucPosition().y() << " "
00660 << (*iter_ExtTrack)->mucPosition().z() << " r "
00661 << (*iter_ExtTrack)->mucPosition().r() << endreq;
00662
00663 if((*iter_ExtTrack)->mucPosition().x() == -99 &&
00664 (*iter_ExtTrack)->mucPosition().y() == -99 &&
00665 (*iter_ExtTrack)->mucPosition().z() == -99)
00666 {
00667 log << MSG::INFO <<"Bad ExtTrack, trackIndex: " << trackIndex << ", skip" << endreq;
00668 continue;
00669 }
00670
00671
00672 krechi = (*iter_ExtTrack)->MucKalchi2();
00673 kdof= (*iter_ExtTrack)->MucKaldof();
00674 kdep = (*iter_ExtTrack)->MucKaldepth();
00675 kbrLay = (*iter_ExtTrack)->MucKalbrLastLayer();
00676 kecLay = (*iter_ExtTrack)->MucKalecLastLayer();
00677 if(kdof<=0)krechi = 0.;
00678 else krechi = krechi/kdof;
00679 kdep = kdep/10.;
00680
00681 RecMucTrack *aTrack = new RecMucTrack();
00682 aTrack->setTrackId(trackIndex);
00683 aTrack->setId(muctrackId);
00684
00685 aTrack->setType(0);
00686
00687 aTrack->SetExtTrack(*iter_ExtTrack);
00688
00689
00690 aTrack->setkalRechi2(krechi);
00691 aTrack->setkalDof(kdof);
00692 aTrack->setkalDepth(kdep);
00693 aTrack->setkalbrLastLayer(kbrLay);
00694 aTrack->setkalecLastLayer(kecLay);
00695
00696
00697 Hep3Vector mdcMomentum = (*iter_MdcTrack)->p3();
00698 Hep3Vector mucPos = (*iter_ExtTrack)->mucPosition();
00699 Hep3Vector mucMomentum = (*iter_ExtTrack)->mucMomentum();
00700
00701
00702 aTrack->SetMdcMomentum( mdcMomentum.x(), mdcMomentum.y(), mdcMomentum.z() );
00703 aTrack->SetExtMucPos( mucPos.x(), mucPos.y(), mucPos.z() );
00704
00705 aTrack->SetExtMucMomentum( mucMomentum.x(), mucMomentum.y(), mucMomentum.z() );
00706 aTrack->SetMucPos( mucPos.x(), mucPos.y(), mucPos.z() );
00707 aTrack->SetMucMomentum( mucMomentum.x(), mucMomentum.y(), mucMomentum.z() );
00708 aTrack->SetCurrentPos( mucPos.x(), mucPos.y(), mucPos.z());
00709 aTrack->SetCurrentDir( mucMomentum.x(), mucMomentum.y(), mucMomentum.z() );
00710 aTrack->SetRecMode(0);
00711
00712 log << MSG::INFO << "Add track, trkIndex: " << trackIndex << "\tmucTrkId: " << muctrackId << endreq;
00713 aRecMucTrackCol->add(aTrack);
00714 muctrackId ++ ;
00715 }
00716 }
00717 else if(m_ExtTrackSeedMode == 3)
00718 {
00719
00720 if (!aMdcTrackCol) {
00721 log << MSG::WARNING << "Can't find MdcTrackCol in TDS!" << endreq;
00722 return StatusCode::SUCCESS;
00723
00724 }
00725 log<< MSG::INFO << "Mdc track size: "<<aMdcTrackCol->size()<<endreq;
00726
00727 int trackIndex = -99;
00728 for(RecMdcTrackCol::iterator iter_mdc1 = aMdcTrackCol->begin(); iter_mdc1 != aMdcTrackCol->end(); iter_mdc1++){
00729
00730
00731
00732
00733 trackIndex = (*iter_mdc1)->trackId();
00734 HepVector helix = (*iter_mdc1)->helix();
00735
00736
00737 float x0 = helix[0] * cos(helix[1]);
00738 float y0 = helix[0] * sin(helix[1]);
00739 float z0 = helix[3];
00740
00741
00742
00743
00744
00745 float dx = -1* sin(helix[1]);
00746 float dy = cos(helix[1]);
00747 float dz = helix[4];
00748
00749
00750
00751 RecMucTrack *aTrack = new RecMucTrack();
00752 aTrack->setTrackId(trackIndex);
00753 aTrack->setId(muctrackId);
00754 muctrackId ++ ;
00755
00756 aTrack->setType(3);
00757
00758
00759 Hep3Vector mucPos(x0,y0,z0);
00760 Hep3Vector mucMomentum(dx,dy,dz);
00761
00762 aTrack->SetExtMucPos(x0,y0,z0);
00763 aTrack->SetExtMucMomentum(dx,dy,dz);
00764
00765 aTrack->SetMucPos(x0,y0,z0);
00766 aTrack->SetMucMomentum(dx,dy,dz);
00767 aTrack->SetCurrentPos(x0,y0,z0);
00768 aTrack->SetCurrentDir(dx,dy,dz);
00769 aTrack->SetRecMode(0);
00770
00771 log << MSG::INFO << "Add track, trkIndex: " << trackIndex << "\tmucTrkId: " << muctrackId << endreq;
00772 aRecMucTrackCol->add(aTrack);
00773 muctrackId ++ ;
00774 }
00775 }
00776 else {log << MSG::INFO <<"ExtTrackSeedMode error!"<<endreq; }
00777
00778
00779 if(m_EmcShowerSeedMode == 1)
00780 {
00781 int trackIndex = 999;
00782 RecEmcShowerCol::iterator iShowerCol;
00783 for(iShowerCol=aShowerCol->begin(); iShowerCol!=aShowerCol->end(); iShowerCol++)
00784 {
00785
00786
00787
00788
00789 RecMucTrack *aTrack = new RecMucTrack();
00790 aTrack->setTrackId(trackIndex);
00791 aTrack->setId(muctrackId);
00792 aTrack->setType(1);
00793
00794 Hep3Vector mucPos = (*iShowerCol)->position();
00795 Hep3Vector mucMomentum = (*iShowerCol)->position();
00796
00797 aTrack->SetExtMucPos( mucPos.x(), mucPos.y(), mucPos.z() );
00798
00799 aTrack->SetExtMucMomentum( mucMomentum.x(), mucMomentum.y(), mucMomentum.z() );
00800 aTrack->SetMucPos( mucPos.x(), mucPos.y(), mucPos.z() );
00801 aTrack->SetMucMomentum( mucMomentum.x(), mucMomentum.y(), mucMomentum.z() );
00802 aTrack->SetCurrentPos( mucPos.x(), mucPos.y(), mucPos.z());
00803 aTrack->SetCurrentDir( mucMomentum.x(), mucMomentum.y(), mucMomentum.z() );
00804 aTrack->SetRecMode(1);
00805
00806
00807 log << MSG::INFO << "Add track, trkIndex: " << trackIndex << "\tmucTrkId: " << muctrackId << endreq;
00808 aRecMucTrackCol->add(aTrack);
00809 muctrackId ++ ;
00810
00811 m_emcrec = 1;
00812 }
00813 }
00814 log << MSG::DEBUG << " track container filled " << endreq;
00815
00816
00817 log << MSG::INFO << "Start tracking..." << endreq;
00818 int runVerbose = 1;
00819 RecMucTrack *aTrack = 0;
00820 for(int iTrack = 0; iTrack < (int)aRecMucTrackCol->size(); iTrack++)
00821 {
00822 log << MSG::DEBUG << "iTrack " << iTrack << endreq;
00823 aTrack = (*aRecMucTrackCol)[iTrack];
00824
00825
00826 Hep3Vector currentPos = aTrack->GetCurrentPos();
00827 Hep3Vector currentDir = aTrack->GetCurrentDir();
00828 if(currentPos.mag() < kMinor) {
00829 log << MSG::WARNING << "No MUC intersection in track " << iTrack << endreq;
00830 continue;
00831 }
00832
00833
00834 int firstHitFound[2] = { 0, 0};
00835 int firstHitGap[2] = {-1, -1};
00836 for(int partSeq = 0; partSeq < (int)MucID::getPartNum(); partSeq++)
00837 {
00838 int iPart = kPartSeq[partSeq];
00839 for(int iGap = 0; iGap < (int)MucID::getGapNum(iPart); iGap++)
00840 {
00841
00842 int seg = -1;
00843 int orient = MucGeoGeneral::Instance()->GetGap(iPart, 0, iGap)->Orient();;
00844
00845 float xInsct, yInsct, zInsct;
00846 aTrack->Project(iPart, iGap, xInsct, yInsct, zInsct, seg);
00847 if(m_MsOutput) cout << "part " << iPart << " gap " << iGap << " x " << xInsct << " y " << yInsct << " z " << zInsct << " seg " << seg << endl;
00848
00849 if(seg == -1) continue;
00850
00851 aTrack->SetCurrentInsct(xInsct, yInsct, zInsct);
00852
00853 for(int iDeltaSeg = 0; iDeltaSeg < kNSegSearch; iDeltaSeg++)
00854 {
00855 int iSeg = seg + kDeltaSeg[iDeltaSeg];
00856 if(iSeg < 0) iSeg += MucID::getSegNum(iPart);
00857 if(iSeg > (int)MucID::getSegNum(iPart) - 1) iSeg -= MucID::getSegNum(iPart);
00858
00859
00860
00861
00862
00863
00864 Hep3Vector mom_mdc = aTrack->getMdcMomentum();
00865 float window = getWindowSize(mom_mdc, iPart, iSeg, iGap);
00866
00867 if (firstHitFound[orient] != 1) window *= kFirstHitWindowRatio;
00868 for (int iHit = 0; iHit < m_MucRecHitContainer->GetGapHitCount(iPart, iSeg, iGap); iHit++)
00869 {
00870 log << MSG::DEBUG << "iSeg " << iSeg << " iHit " << iHit << endreq;
00871 MucRecHit* pHit = m_MucRecHitContainer->GetHit(iPart, iSeg, iGap, iHit);
00872
00873
00874 if (!pHit) {
00875 log << MSG::WARNING << "MucRecTrkExt: null pointer to pHit" << endreq;
00876 continue;
00877 }
00878 else
00879 {
00880
00881
00882
00883 float dX = aTrack->GetHitDistance2(pHit);
00884 log << MSG::DEBUG << "distance = " << setw(8) << dX << " size " << setw(4) << window << endreq;
00885
00886 if(m_NtOutput >= 2){
00887 m_distance = dX;
00888 m_part = iPart; m_seg = iSeg; m_gap = iGap; m_strip = pHit->Strip();
00889 m_diff = -99; m_tuple->write();
00890 }
00891
00892 if (fabs(dX) < window)
00893 {
00894
00895
00896
00897
00898 if(aTrack->GetRecMode() == 0) {
00899 pHit->SetHitMode(1);
00900 aTrack->AttachHit(pHit);
00901
00902 }
00903 if(aTrack->GetRecMode() == 1) {
00904
00905 if(pHit->GetHitMode()!=1) {
00906 aTrack->AttachHit(pHit);
00907 pHit->SetHitMode(2);
00908 }
00909 }
00910
00911
00912 aTrack->pushExtDistHits(dX);
00913
00914 if (firstHitGap[orient] == -1) firstHitGap[orient] = iGap;
00915 firstHitFound[orient] = 1;
00916
00917
00918 log << MSG::DEBUG << " part " << iPart << " seg " << iSeg << " gap " << iGap
00919 << " strip " << setw(2) << pHit->Strip() << " attatched" << endreq;
00920 log << MSG::DEBUG << "current total hits " << aTrack->GetTotalHits() << endreq;
00921 }
00922 else {
00923 m_NHitsLostInGap[iGap]++;
00924
00925
00926
00927 }
00928 }
00929 }
00930 aTrack->CalculateInsct(iPart, iSeg, iGap);
00931 }
00932
00933
00934 if(m_ExtTrackSeedMode != 3 && !m_Blind) {
00935 if (firstHitFound[orient] && firstHitGap[orient] != iGap) aTrack->CorrectDir();
00936 aTrack->CorrectPos();
00937 }
00938
00939
00940 aTrack->AttachInsct(aTrack->GetCurrentInsct());
00941 }
00942 }
00943 aTrack->LineFit(m_fittingMethod);
00944 aTrack->ComputeTrackInfo(m_fittingMethod);
00945 log << MSG::INFO <<"Fit track done! trackIndex: " << aTrack->trackId() << ", mucId: " << aTrack->id() << ", RecMode: " << aTrack->GetRecMode() << endreq;
00946
00947 if(m_NtOutput>=1)
00948 {
00949 m_depth = aTrack->depth();
00950 m_Distance = aTrack->distance();
00951 m_MaxHits = aTrack->maxHitsInLayer();
00952 m_Chi2 = aTrack->chi2();
00953 m_Dist_x = aTrack->GetExtMucDistance().x();
00954 m_Dist_y = aTrack->GetExtMucDistance().y();
00955 m_Dist_z = aTrack->GetExtMucDistance().z();
00956 m_posx_ext = aTrack->GetMucStripPos().x();
00957 m_posy_ext = aTrack->GetMucStripPos().y();
00958 m_posz_ext = aTrack->GetMucStripPos().z();
00959
00960 m_emctrack = m_emcrec;
00961 }
00962
00963
00964 if(m_NtOutput>=2)
00965 {
00966 vector<MucRecHit*> attachedHits = aTrack->GetHits();
00967 vector<MucRecHit*> expectedHits = aTrack->GetExpectedHits();
00968 vector<float> distanceHits = aTrack->getDistHits();
00969 vector<float> distanceHits_quad = aTrack->getQuadDistHits();
00970 vector<float> distanceHits_ext = aTrack->getExtDistHits();
00971
00972 for(int i=0; i< expectedHits.size(); i++)
00973 {
00974 MucRecHit *ihit = expectedHits[i];
00975
00976 }
00977
00978 for(int j=0; j< attachedHits.size(); j++)
00979 {
00980 MucRecHit *jhit = attachedHits[j];
00981
00982 if(attachedHits.size()==distanceHits.size())
00983 {
00984 m_part = jhit->Part(); m_seg = jhit->Seg(); m_gap = jhit->Gap();
00985 m_strip = jhit->Strip();
00986 m_distance = distanceHits[j];
00987 m_Distance = distanceHits_quad[j];
00988 m_diff = -9999;
00989
00990
00991
00992 m_tuple->write();
00993 }
00994 }
00995 }
00996
00997 int nHitsAttached = aTrack->GetTotalHits();
00998
00999
01000 log << MSG::DEBUG << "track Index " << aTrack->trackId()
01001 << setw(2) << mucDigiCol->size() - nHitsAttached << " of "
01002 << setw(2) << mucDigiCol->size() << " lost " << endreq;
01003
01004 }
01005
01006
01007
01008
01009
01010
01011
01012 if(m_MucHitSeedMode == 1)
01013 {
01014 MucRecHit *pHit = 0 ,*pHit0 = 0, *pHit1 = 0;
01015 int count, orient;
01016
01017 for (int iPart = 0; iPart < (int)MucID::getPartNum(); iPart++)
01018 {
01019 for (int iSeg = 0; iSeg < (int)MucID::getSegNum(iPart); iSeg++)
01020 {
01021 bool hit0 = false, hit1 = false; int firstgap0 = -1, firstgap1 = -1; int nStrip0 = 0, nStrip1 = 0;
01022 Hep3Vector posHit0, posHit1;
01023 for(int iGap = 0; iGap < (int)MucID::getGapNum(iPart); iGap++)
01024 {
01025 count = m_MucRecHitContainer->GetGapHitCount(iPart, iSeg, iGap);
01026 for (int iHit0 = 0; iHit0 < count; iHit0++)
01027 {
01028
01029 pHit = m_MucRecHitContainer->GetHit(iPart, iSeg, iGap, iHit0);
01030 if (!pHit) {
01031 log << MSG::FATAL << "MucRecRoadFinder-E10 " << " part " << iPart
01032 << " seg " << iSeg << " gap " << iGap << " null pointer"
01033 << endl;
01034 }
01035 else
01036 {
01037
01038 if(pHit->GetHitMode() == -1)
01039 {
01040 orient = MucGeoGeneral::Instance()->GetGap(iPart, 0, iGap)->Orient();
01041 if(orient == 1 && hit0 == false){
01042 hit0 = true;
01043 firstgap0 = iGap;
01044 }
01045 if(iGap == firstgap0){
01046 posHit0 += pHit->GetCenterPos();
01047 nStrip0++;
01048 }
01049
01050 if(orient == 0 && hit1 == false){
01051 hit1 = true;
01052 firstgap1 = iGap;
01053 }
01054 if(iGap == firstgap1){
01055 posHit1 += pHit->GetCenterPos();
01056 nStrip1++;
01057 }
01058
01059
01060
01061
01062
01063
01064
01065
01066 }
01067 }
01068 }
01069 }
01070
01071
01072 if(hit0 && hit1)
01073 {
01074 posHit0 /= nStrip0;
01075 posHit1 /= nStrip1;
01076
01077
01078
01079 int trackIndex = 999;
01080 RecMucTrack *aTrack = new RecMucTrack();
01081 aTrack->setTrackId(trackIndex);
01082 aTrack->setId(muctrackId);
01083 muctrackId ++ ;
01084
01085 aTrack->setType(2);
01086
01087 Hep3Vector mucPos, mucMomentum;
01088 if(iPart==1) {
01089 mucMomentum.set(posHit0.x(),posHit0.y(),posHit1.z());
01090 mucPos = mucMomentum * 0.9;
01091 }
01092 else{
01093 mucMomentum.set(posHit0.x(),posHit1.y(),posHit0.z()*0.5+posHit1.z()*0.5);
01094 mucPos = mucMomentum * 0.9;
01095 }
01096
01097 aTrack->SetExtMucPos( mucPos.x(), mucPos.y(), mucPos.z() );
01098
01099 aTrack->SetExtMucMomentum( mucMomentum.x(), mucMomentum.y(), mucMomentum.z() );
01100 aTrack->SetMucPos( mucPos.x(), mucPos.y(), mucPos.z() );
01101 aTrack->SetMucMomentum( mucMomentum.x(), mucMomentum.y(), mucMomentum.z() );
01102 aTrack->SetCurrentPos( mucPos.x(), mucPos.y(), mucPos.z());
01103 aTrack->SetCurrentDir( mucMomentum.x(), mucMomentum.y(), mucMomentum.z() );
01104
01105 aTrack->SetRecMode(2);
01106 aRecMucTrackCol->add(aTrack);
01107
01108 TrackFinding(aTrack);
01109 }
01110
01111 if(m_NtOutput>=1){
01112 m_depth = aTrack->depth();
01113 m_Distance = aTrack->distance();
01114 m_MaxHits = aTrack->maxHitsInLayer();
01115 m_Chi2 = aTrack->chi2();
01116 m_Dist_x = aTrack->GetExtMucDistance().x();
01117 m_Dist_y = aTrack->GetExtMucDistance().y();
01118 m_Dist_z = aTrack->GetExtMucDistance().z();
01119 m_posx_ext = aTrack->GetMucStripPos().x();
01120 m_posy_ext = aTrack->GetMucStripPos().y();
01121 m_posz_ext = aTrack->GetMucStripPos().z();
01122
01123 m_emctrack = 2;
01124 }
01125
01126 }
01127 }
01128 }
01129
01130
01131 int nTracksTotal = 0;
01132 int nTracksFound = 0;
01133 int nTracksLost = 0;
01134 int nTracksLostByExt = 0;
01135 int nTracksMisFound = 0;
01136
01137 int nDigisTotal = 0;
01138 int nHitsTotal = 0;
01139 int nHitsFound = 0;
01140 int nHitsLost = 0;
01141 int nHitsMisFound = 0;
01142
01143
01144
01145
01146
01147
01148
01149
01150
01151
01152
01153
01154
01155
01156
01157
01158
01159
01160
01161
01162
01163
01164
01165
01166
01167
01168
01169
01170
01171
01172
01173
01174
01175
01176
01177
01178
01179
01180
01181
01182
01183
01184
01185
01186
01187
01188
01189
01190
01191
01192
01193
01194
01195
01196
01197
01198
01199
01200
01201
01202
01203
01204
01205
01206
01207
01208
01209
01210
01211
01212
01213
01214
01215
01216
01217
01218
01219
01220
01221
01222
01223
01224
01225
01226
01227
01228
01229
01230
01231
01232
01233
01234
01235
01236
01237
01238
01239
01240
01241
01242
01243
01244
01245
01246
01247
01248
01249
01250
01251
01252
01253
01254
01255
01256
01257
01258
01259
01260
01261
01262
01263
01264
01265
01266
01267
01268
01269
01270
01271
01272
01273
01274
01275
01276
01277
01278
01279
01280
01281
01282
01283
01284
01285
01286
01287
01288
01289
01290
01291
01292
01293
01294
01295
01296
01297
01298
01299
01300
01301
01302
01303
01304
01305
01306
01307
01308
01309
01310
01311
01312
01313
01314
01315
01316
01317
01318
01319
01320
01321
01322
01323
01324
01325
01326
01327
01328
01329
01330
01331
01332
01333
01334
01335
01336
01337
01338
01339
01340
01341
01342
01343
01344
01345
01346
01347 m_NDigisTotal += nDigisTotal;
01348 m_NHitsTotal += nHitsTotal;
01349 m_NHitsFoundTotal += nHitsFound;
01350 m_NHitsLostTotal += nHitsLost;
01351 m_NHitsMisFoundTotal += nHitsMisFound;
01352
01353
01354 m_NTracksTotal += nTracksTotal;
01355 m_NTracksRecoTotal += nTracksFound;
01356 m_NTracksLostByMdcTotal += nTracksLost;
01357 m_NTracksLostByExtTotal += nTracksLostByExt;
01358 m_NTracksMdcGhostTotal += nTracksMisFound;
01359 if (aRecMucTrackCol->size() > 0)
01360 {
01361 RecMucTrack *aRecMucTrack = (*aRecMucTrackCol)[0];
01362
01363 if(m_NtOutput>=1){
01364 m_posx = aRecMucTrack->getMucPos().x();
01365 m_posy = aRecMucTrack->getMucPos().y();
01366 m_posz = aRecMucTrack->getMucPos().z();
01367
01368 m_posx_sigma = aRecMucTrack->getMucPosSigma().x();
01369 m_posy_sigma = aRecMucTrack->getMucPosSigma().y();
01370 m_posz_sigma = aRecMucTrack->getMucPosSigma().z();
01371 }
01372
01373
01374 }
01375
01376 if(m_NtOutput>=1) m_tuple->write();
01377 return StatusCode::SUCCESS;
01378 }
01379
01380
01381 StatusCode MucRecTrkExt::finalize() {
01382
01383 MsgStream log(msgSvc(), name());
01384 log << MSG::INFO << "in finalize()" << endreq;
01385
01386 log << MSG::INFO << " In " << m_NDigisTotal << " total MucDigis " << endreq;
01387 log << MSG::INFO << " In " << m_NHitsTotal << " total MucMcHits " << endreq;
01388 log << MSG::INFO << m_NHitsFoundTotal << " hits found in Reco tracks, " << endreq;
01389 log << MSG::INFO << m_NHitsLostTotal << " hits lost (not found), of which "
01390 << m_NHitsLostByMdcTotal << " lost because Mdc track not constructed "
01391 << m_NHitsLostByExtTotal << " lost because Ext track not intersect with muc " << endreq;
01392 log << MSG::INFO << m_NHitsMisFoundTotal << " hits mis found (not belong to this track, but mis attached)" << endreq;
01393
01394 log << MSG::INFO << " In " << m_NTracksTotal << " total Mc tracks " << endreq;
01395 log << MSG::INFO << m_NTracksRecoTotal << " tracks found " << endreq;
01396 log << MSG::INFO << m_NTracksLostByMdcTotal << " tracks lost because Mdc track not constructed " << endreq;
01397 log << MSG::INFO << m_NTracksLostByExtTotal << " tracks lost because Ext track not intersect with muc " << endreq;
01398 log << MSG::INFO << m_NTracksMdcGhostTotal << " tracks are Mdc ghost tracks " << endreq;
01399
01400
01401
01402
01403
01404
01405
01406
01407
01408
01409
01410
01411
01412
01413
01414
01415
01416
01417 return StatusCode::SUCCESS;
01418 }
01419
01420 void
01421 MucRecTrkExt::TrackFinding(RecMucTrack *aTrack)
01422 {
01423 MsgStream log(msgSvc(), name());
01424
01425 Hep3Vector currentPos = aTrack->GetCurrentPos();
01426 Hep3Vector currentDir = aTrack->GetCurrentDir();
01427
01428
01429
01430
01431
01432 int firstHitFound[2] = { 0, 0};
01433 int firstHitGap[2] = {-1, -1};
01434 for(int partSeq = 0; partSeq < (int)MucID::getPartNum(); partSeq++)
01435 {
01436 int iPart = kPartSeq[partSeq];
01437 for(int iGap = 0; iGap < (int)MucID::getGapNum(iPart); iGap++)
01438 {
01439 int seg = -1;
01440 int orient = MucGeoGeneral::Instance()->GetGap(iPart, 0, iGap)->Orient();;
01441 float xInsct, yInsct, zInsct;
01442 aTrack->Project(iPart, iGap, xInsct, yInsct, zInsct, seg);
01443 log << MSG::INFO << "part " << iPart << " gap " << iGap << " x " << xInsct << " y " << yInsct << " z " << zInsct << " seg " << seg << endreq;
01444
01445
01446 if(seg == -1) {
01447
01448
01449 continue;
01450 }
01451 aTrack->SetCurrentInsct(xInsct, yInsct, zInsct);
01452
01453 for(int iDeltaSeg = 0; iDeltaSeg < kNSegSearch; iDeltaSeg++)
01454 {
01455 int iSeg = seg + kDeltaSeg[iDeltaSeg];
01456 if(iSeg < 0) iSeg += MucID::getSegNum(iPart);
01457 if(iSeg > (int)MucID::getSegNum(iPart) - 1) iSeg -= MucID::getSegNum(iPart);
01458
01459
01460 Hep3Vector mom_mdc = aTrack->getMdcMomentum();
01461 float window = getWindowSize(mom_mdc, iPart, iSeg, iGap);
01462
01463 if (firstHitFound[orient] != 1) window *= kFirstHitWindowRatio;
01464
01465 for (int iHit = 0; iHit < m_MucRecHitContainer->GetGapHitCount(iPart, iSeg, iGap); iHit++)
01466 {
01467 log << MSG::DEBUG << "iSeg " << iSeg << " iHit " << iHit << endreq;
01468 MucRecHit* pHit = m_MucRecHitContainer->GetHit(iPart, iSeg, iGap, iHit);
01469
01470
01471 if (!pHit) {
01472 log << MSG::WARNING << "MucRecTrkExt: null pointer to pHit" << endreq;
01473 continue;
01474 }
01475 else
01476 {
01477
01478 float dX = aTrack->GetHitDistance2(pHit);
01479 log << MSG::DEBUG << "distance = " << setw(8) << dX << " size " << setw(4) << window << endreq;
01480
01481
01482 if (fabs(dX) < window)
01483 {
01484
01485
01486
01487
01488 if(aTrack->GetRecMode() == 0) {
01489 pHit->SetHitMode(1);
01490 aTrack->AttachHit(pHit);
01491
01492 }
01493 if(aTrack->GetRecMode() == 1) {
01494
01495 if(pHit->GetHitMode()!=1) {
01496 aTrack->AttachHit(pHit);
01497 pHit->SetHitMode(2);
01498 }
01499 }
01500 if(aTrack->GetRecMode() == 2) {
01501
01502 if(pHit->GetHitMode()==-1) {
01503 aTrack->AttachHit(pHit);
01504 pHit->SetHitMode(2);
01505 }
01506 }
01507
01508 if (firstHitGap[orient] == -1) firstHitGap[orient] = iGap;
01509 firstHitFound[orient] = 1;
01510
01511
01512
01513 log << MSG::DEBUG << " part " << iPart << " seg " << iSeg << " gap " << iGap
01514 << " strip " << setw(2) << pHit->Strip() << " attatched" << endreq;
01515 log << MSG::DEBUG << "current total hits " << aTrack->GetTotalHits() << endreq;
01516 }
01517 else
01518 {
01519 m_NHitsLostInGap[iGap]++;
01520
01521
01522
01523 }
01524 }
01525 }
01526
01527 aTrack->CalculateInsct(iPart, iSeg, iGap);
01528 }
01529
01530
01531 if (firstHitFound[orient] && firstHitGap[orient] != iGap) aTrack->CorrectDir();
01532 aTrack->CorrectPos();
01533
01534
01535 aTrack->AttachInsct(aTrack->GetCurrentInsct());
01536 }
01537 }
01538
01539 aTrack->LineFit(m_fittingMethod);
01540 aTrack->ComputeTrackInfo(m_fittingMethod);
01541
01542
01543 }
01544
01545
01546 float
01547 MucRecTrkExt::getWindowSize(Hep3Vector mom, int part, int seg, int gap)
01548 {
01549 int mom_id;
01550 if(mom.mag()<0.6) mom_id = 0;
01551 else if(mom.mag()<0.8) mom_id = 1;
01552 else if(mom.mag()<1) mom_id = 2;
01553 else mom_id = 3;
01554
01555 return kWindowSize_Mom_Gap[mom_id][gap];
01556 }
01557
01558