Definition at line 171 of file MucRecRoadFinder.cxx.
References MucRecHitContainer::AddHit(), aMucRecHitContainer, MucRec2DRoad::AttachHit(), MucRec2DRoad::AttachHitNoFit(), MucID::barrel_ec(), MucID::channel(), MucRecHitContainer::Clear(), count, Bes_Common::DEBUG, calibUtil::ERROR, Bes_Common::FATAL, MucRecHit::Gap(), RecMucTrack::getDistHits(), RecMucTrack::GetExpectedHits(), MucRecHit::GetGap(), MucRecHitContainer::GetGapHitCount(), MucRecHitContainer::GetHit(), MucRec2DRoad::GetHitDistance(), MucRecHit::GetHitMode(), RecMucTrack::GetHits(), MucRec2DRoad::GetHits(), MucRec3DRoad::GetLastGap(), MucRec2DRoad::GetLastGap(), MucRec3DRoad::GetLastGapDelta(), RecMucTrack::getMucMomentum(), RecMucTrack::getMucPos(), MucRecHitContainer::GetMucRecHitCol(), MucRec2DRoad::GetNGapsWithHits(), MucRec3DRoad::GetPart(), MucRec2DRoad::GetPart(), MucID::getPartNum(), MucRec3DRoad::GetSeg(), MucID::getSegNum(), MucRec2DRoad::GetSlope(), MucRec3DRoad::GetTotalHits(), MucRec2DRoad::GetTotalHits(), MucRec3DRoad::GetTotalHitsDelta(), genRecEmupikp::i, Bes_Common::INFO, ganga-rec::j, kMaxDeltaLastGap, kMaxDeltaTotalHits, kMaxSkippedGaps, kMaxXChisq, kMaxYChisq, kMinFiredGaps, kMinLastGap2D, kMinLastGap3D, kMinSharedHits2D, kNSeedLoops, kSearchLength, kSearchOrder, kWindowSize, MucID::layer(), m_angle_cosmic, m_angle_updown, m_diff, m_dist, m_emcDown, m_emcUp, m_event, m_filter_event, m_gap, m_maxhit, m_maxHitsRec, m_mccosmic, m_MsOutput, m_mucDown, m_mucUp, m_multihit, m_NEvent, m_NEventReco, m_ngapwithhits, m_nhit, m_NtOutput, m_onlyseedfit, m_part, m_phi, m_phi_mc, m_phi_mc_pipe, m_phi_pipe, m_projx, m_projz, m_px, m_px_mc, m_py, m_py_mc, m_pz, m_pz_mc, m_run, m_seedtype, m_seg, m_strip, m_theta, m_theta_mc, m_theta_mc_pipe, m_theta_pipe, m_tuple, m_united, DstMucTrack::maxHitsInLayer(), msgSvc(), EventModel::Recon::MucRecHitCol, DstMucTrack::numHits(), DstMucTrack::numLayers(), MucGeoGap::Orient(), MucRecHit::Part(), pid, MucRec3DRoad::ProjectWithSigma(), EventModel::Recon::RecMucTrackCol, release, MucRecHit::Seg(), MucID::segment(), RecMucTrack::SetCurrentDir(), RecMucTrack::SetCurrentPos(), RecMucTrack::SetExtMucMomentum(), RecMucTrack::SetExtMucPos(), MucRecHit::SetHitMode(), MucRecHit::SetHitSeed(), DstMucTrack::setId(), MucRec2DRoad::SetMaxNSkippedGaps(), RecMucTrack::SetMucMomentum(), RecMucTrack::SetMucPos(), MucRecHitContainer::SetMucRecHitCol(), RecMucTrack::SetRecMode(), MucRec3DRoad::SetSimpleFitParams(), RecMucTrack::setTrackId(), sign(), deljobs::string, MucRecHit::Strip(), TrackFinding(), DstMucTrack::trackId(), MucRec3DRoad::TransformPhiRToXY(), Bes_Common::WARNING, and x.
00171 {
00172
00173 MsgStream log(msgSvc(), name());
00174 log << MSG::INFO << "in execute()" << endreq;
00175
00176
00177 int event, run;
00178
00179 SmartDataPtr<Event::EventHeader> eventHeader(eventSvc(),"/Event/EventHeader");
00180 if (!eventHeader) {
00181 log << MSG::FATAL << "Could not find Event Header" << endreq;
00182 return( StatusCode::FAILURE);
00183 }
00184 log << MSG::INFO << "Run: " << eventHeader->runNumber() << " Event: " << eventHeader->eventNumber() << endreq;
00185
00186 event = eventHeader->eventNumber();
00187 run = eventHeader->runNumber();
00188
00189
00190 string release = getenv("BES_RELEASE");
00191
00192
00193
00194 for (std::vector<FilterEvent>::iterator it = m_filter_event.begin();
00195 it != m_filter_event.end(); ++it) {
00196 const FilterEvent& fe = (*it);
00197 if (release == fe.bossver && run == fe.runid && event == fe.eventid) {
00198 cout << "SKIP: " << fe.bossver << " "
00199 << fe.runid << " "
00200 << fe.eventid << std::endl;
00201 return StatusCode::SUCCESS;
00202 }
00203 }
00204
00205 int digiId;
00206
00207
00208
00209 Hep3Vector cosmicMom;
00210 if(m_mccosmic==1){
00211 SmartDataPtr<McParticleCol> mcParticleCol(eventSvc(),"/Event/MC/McParticleCol");
00212 if (!mcParticleCol) {
00213 log << MSG::FATAL << "Could not find McParticle" << endreq;
00214 return( StatusCode::FAILURE);
00215 }
00216
00217 McParticleCol::iterator iter_mc = mcParticleCol->begin();
00218 for (;iter_mc != mcParticleCol->end(); iter_mc++, digiId++) {
00219 log << MSG::DEBUG
00220 << " particleId = " << (*iter_mc)->particleProperty()
00221 << endreq;
00222 int pid = (*iter_mc)->particleProperty();
00223
00224 if(fabs(pid)!=13) continue;
00225
00226 HepLorentzVector initialMomentum = (*iter_mc)->initialFourMomentum();
00227 Hep3Vector startP(initialMomentum.px(), initialMomentum.py(), initialMomentum.pz());
00228 Hep3Vector rotate(-initialMomentum.px(),initialMomentum.pz(),initialMomentum.py());
00229
00230 if(m_NtOutput>=1){
00231 m_px_mc = initialMomentum.px();
00232 m_py_mc = initialMomentum.py();
00233 m_pz_mc = initialMomentum.pz();
00234
00235 m_theta_mc = rotate.theta();
00236 m_phi_mc = rotate.phi();
00237
00238 m_theta_mc_pipe = startP.theta();
00239 m_phi_mc_pipe = startP.phi();
00240
00241
00242 }
00243
00244
00245 cosmicMom = startP;
00246
00247 }
00248 }
00249
00250
00251
00252
00253
00254
00255
00256
00257
00258
00259
00260
00261
00262
00263
00264
00265
00266
00267
00268
00269
00270
00271
00272
00273
00274
00275
00276
00277
00278
00279
00280
00281
00282
00283
00284
00285
00286
00287
00288
00289
00290
00291
00292
00293
00294
00295
00296
00297
00298
00299
00300
00301
00302
00303
00304
00305 SmartDataPtr<MucDigiCol> mucDigiCol(eventSvc(),"/Event/Digi/MucDigiCol");
00306 if (!mucDigiCol) {
00307 log << MSG::FATAL << "Could not find MUC digi" << endreq;
00308 return( StatusCode::FAILURE);
00309 }
00310
00311 MucDigiCol::iterator iter4 = mucDigiCol->begin();
00312 digiId = 0;
00313 for (;iter4 != mucDigiCol->end(); iter4++, digiId++) {
00314 }
00315 log << MSG::INFO << "MUC digis:: " << digiId << endreq;
00316 if( digiId == 0) return( StatusCode::SUCCESS );
00317
00318
00319 RecMucTrackCol *aMucTrackCol;
00320 int trackId = -1;
00321 int muctrackId = -1;
00322
00323 if(m_united != 1)
00324 {
00325 Identifier mucID;
00326 if (!aMucRecHitContainer) {
00327 aMucRecHitContainer = new MucRecHitContainer();
00328 }
00329 aMucRecHitContainer->Clear();
00330
00331 MucRecHitCol *aMucRecHitCol = new MucRecHitCol();
00332 aMucRecHitContainer->SetMucRecHitCol(aMucRecHitCol);
00333
00334
00335 SmartIF<IDataManagerSvc> dataManSvc(eventSvc());
00336 DataObject *mucRecHitCol;
00337 eventSvc()->findObject("/Event/Recon/MucRecHitCol",mucRecHitCol);
00338 if(mucRecHitCol != NULL) {
00339 dataManSvc->clearSubTree("/Event/Recon/MucRecHitCol");
00340 eventSvc()->unregisterObject("/Event/Recon/MucRecHitCol");
00341 }
00342
00343 StatusCode sc = eventSvc()->registerObject("/Event/Recon/MucRecHitCol", aMucRecHitContainer->GetMucRecHitCol());
00344
00345 MucDigiCol::iterator iterMuc = mucDigiCol->begin();
00346 int mucDigiId = 0;
00347 for (;iterMuc != mucDigiCol->end(); iterMuc++, mucDigiId++) {
00348 mucID = (*iterMuc)->identify();
00349 int part = MucID::barrel_ec(mucID);
00350 int seg = MucID::segment(mucID);
00351 int gap = MucID::layer(mucID);
00352 int strip = MucID::channel(mucID);
00353
00354
00355
00356 aMucRecHitContainer->AddHit(part, seg, gap, strip);
00357 log << MSG::INFO << " digit" << mucDigiId << " : "
00358 << " part " << part
00359 << " seg " << seg
00360 << " gap " << gap
00361 << " strip " << strip
00362 << endreq;
00363 }
00364
00365 DataObject *aReconEvent ;
00366 eventSvc()->findObject("/Event/Recon",aReconEvent);
00367 if(aReconEvent==NULL) {
00368
00369 aReconEvent = new ReconEvent();
00370 StatusCode sc = eventSvc()->registerObject("/Event/Recon",aReconEvent);
00371 if(sc!=StatusCode::SUCCESS) {
00372 log << MSG::FATAL << "Could not register ReconEvent" <<endreq;
00373 return( StatusCode::FAILURE);
00374 }
00375 }
00376 StatusCode fr = eventSvc()->findObject("/Event/Recon",aReconEvent);
00377 if(fr!=StatusCode::SUCCESS) {
00378 log << MSG::WARNING << "Could not find register ReconEvent, will create it" <<endreq;
00379 StatusCode sc = eventSvc()->registerObject("/Event/Recon",aReconEvent);
00380 if(sc!=StatusCode::SUCCESS) {
00381 log << MSG::FATAL << "Could not register ReconEvent" <<endreq;
00382 return( StatusCode::FAILURE);
00383 }
00384 }
00385
00386
00387
00388 DataObject *mucTrackCol;
00389 eventSvc()->findObject("/Event/Recon/RecMucTrackCol",mucTrackCol);
00390 if(mucTrackCol != NULL) {
00391 dataManSvc->clearSubTree("/Event/Recon/RecMucTrackCol");
00392 eventSvc()->unregisterObject("/Event/Recon/RecMucTrackCol");
00393 }
00394
00395 aMucTrackCol = new RecMucTrackCol();
00396 sc = eventSvc()->registerObject("/Event/Recon/RecMucTrackCol", aMucTrackCol);
00397 aMucTrackCol->clear();
00398
00399
00400 SmartDataPtr<RecMucTrackCol> findMucTrackCol(eventSvc(),"/Event/Recon/RecMucTrackCol");
00401 if (!findMucTrackCol) {
00402 log << MSG::FATAL << "Could not find RecMucTrackCol" << endreq;
00403 return( StatusCode::FAILURE);
00404 }
00405 aMucTrackCol->clear();
00406
00407
00408 log << MSG::INFO <<"MaxHitsRec : "<<m_maxHitsRec<< endreq;
00409 if(digiId> m_maxHitsRec) return StatusCode::SUCCESS;
00410
00411
00412
00413 trackId = 0;
00414 muctrackId = 0;
00415 }
00416 else if(m_united == 1){
00417
00418 Identifier mucID;
00419 if (!aMucRecHitContainer) {
00420 aMucRecHitContainer = new MucRecHitContainer();
00421 }
00422 aMucRecHitContainer->Clear();
00423
00424 SmartDataPtr<MucRecHitCol> aMucRecHitCol (eventSvc(),"/Event/Recon/MucRecHitCol");
00425 if(aMucRecHitCol == NULL) {
00426 log << MSG::WARNING << "MucRecHitCol is NULL" << endreq;
00427 }
00428
00429 aMucRecHitContainer->SetMucRecHitCol(aMucRecHitCol);
00430
00431 SmartDataPtr<RecMucTrackCol> mucTrackCol (eventSvc(),"/Event/Recon/RecMucTrackCol");
00432 if(mucTrackCol == NULL) {
00433 log << MSG::WARNING << "aMucTrackCol is NULL" << endreq;
00434 }
00435
00436 log << MSG::INFO <<"mucTrackCol size: "<<mucTrackCol->size()<<" hitCol size: "<<aMucRecHitCol->size()<<endreq;
00437 aMucTrackCol = mucTrackCol;
00438
00439
00440 SmartDataPtr<RecExtTrackCol> aExtTrackCol(eventSvc(),"/Event/Recon/RecExtTrackCol");
00441 if (!aExtTrackCol) {
00442 log << MSG::WARNING << "Can't find ExtTrackCol in TDS!" << endreq;
00443 }
00444
00445 SmartDataPtr<RecMdcTrackCol> aMdcTrackCol(eventSvc(),"/Event/Recon/RecMdcTrackCol");
00446 if (!aMdcTrackCol) {
00447 log << MSG::WARNING << "Can't find MdcTrackCol in TDS!" << endreq;
00448 }
00449
00450 if( aExtTrackCol && aMdcTrackCol ) trackId = aMdcTrackCol->size();
00451 muctrackId = aMucTrackCol->size();
00452 }
00453
00454 int hasEmcUp = 0;
00455 int hasEmcDown = 0;
00456 SmartDataPtr<RecEmcShowerCol> aShowerCol(eventSvc(),"/Event/Recon/RecEmcShowerCol");
00457 if (!aShowerCol) {
00458 log << MSG::WARNING << "Could not find RecEmcShowerCol" << endreq;
00459
00460 }
00461 else{
00462 RecEmcShowerCol::iterator iShowerCol;
00463 for(iShowerCol=aShowerCol->begin();
00464 iShowerCol!=aShowerCol->end();
00465 iShowerCol++){
00466
00467
00468
00469
00470
00471
00472 if(((*iShowerCol)->position()).phi()>0&&((*iShowerCol)->position()).phi()<3.1415926) hasEmcUp++;
00473 else hasEmcDown++;
00474 }
00475 }
00476 if(m_NtOutput >= 1){
00477 m_emcUp = hasEmcUp; m_emcDown = hasEmcDown;
00478 }
00479
00480
00481 m_NEvent++;
00482 MucRec3DRoadContainer *p3DRoadC = new MucRec3DRoadContainer();
00483
00484 if (!p3DRoadC) {
00485 p3DRoadC = new MucRec3DRoadContainer();
00486 }
00487 p3DRoadC->clear();
00488
00489 MucRec2DRoadContainer *p2DRoad0C = new MucRec2DRoadContainer();
00490
00491 if (!p2DRoad0C) {
00492 p2DRoad0C = new MucRec2DRoadContainer();
00493 }
00494 p2DRoad0C->clear();
00495
00496 MucRec2DRoadContainer *p2DRoad1C = new MucRec2DRoadContainer();
00497
00498 if (!p2DRoad1C) {
00499 p2DRoad1C = new MucRec2DRoadContainer();
00500 }
00501 p2DRoad1C->clear();
00502 log << MSG::INFO << "Muc 2D 3D Road Container created" << endreq;
00503
00504
00505
00506
00507
00508
00509
00510
00511 MucRecHit *pHit0 = 0, *pHit1 = 0;
00512 int count0, count1, count, iGap0, iGap1;
00513 int orient;
00514
00515 for (int iPart = 0; iPart < (int)MucID::getPartNum(); iPart++) {
00516 for (int iSeg = 0; iSeg < (int)MucID::getSegNum(iPart); iSeg++) {
00517 for (int iOrient = 0; iOrient < 2; iOrient++) {
00518 int nLoops = 1;
00519 if(m_seedtype == 1) nLoops = kNSeedLoops;
00520 for (int iLoop = 0; iLoop < nLoops; iLoop++) {
00521
00522 int length = kSearchLength[iLoop][iPart][iOrient];
00523 if(m_seedtype == 1){
00524 iGap0 = kSearchOrder[iLoop][iPart][iOrient][0];
00525 iGap1 = kSearchOrder[iLoop][iPart][iOrient][1];
00526 }
00527 else {
00528 int setgap0 = 0;
00529 MucDigiCol::iterator iterMuc = mucDigiCol->begin();
00530 int mucDigiId = 0;
00531 Identifier mucID;
00532 iGap0 = 0; iGap1 = 0;
00533 for (;iterMuc != mucDigiCol->end(); iterMuc++, mucDigiId++) {
00534 mucID = (*iterMuc)->identify();
00535 int part = MucID::barrel_ec(mucID);
00536 int seg = MucID::segment(mucID);
00537 int gap = MucID::layer(mucID);
00538 int strip = MucID::channel(mucID);
00539 int orient = 0;
00540 if ( (part == 1 && gap % 2 == 0) || (part != 1 && gap % 2 == 1) ) orient = 1;
00541
00542 if(part == iPart && seg == iSeg && orient == iOrient && setgap0 == 0){
00543 iGap0 = gap;
00544 setgap0 = 1;
00545 }
00546 if(part == iPart && seg == iSeg && orient == iOrient && setgap0 == 1 && gap != iGap0 ){
00547 iGap1 = gap;
00548 }
00549 }
00550 }
00551
00552 if(m_MsOutput) cout <<"Find seed gap in ( "<<iPart<<" "<<iSeg<<" ) gap0: "<<iGap0<<" gap1: "<<iGap1<<endl;
00553
00554
00555 if(iGap0 > iGap1){
00556 int tempGap = iGap0;
00557 iGap0 = iGap1;
00558 iGap1 = tempGap;
00559 }
00560 if(iGap1 == iGap0) continue;
00561
00562 count0 = aMucRecHitContainer->GetGapHitCount(iPart, iSeg, iGap0);
00563 for (int iHit0 = 0; iHit0 < count0; iHit0++) {
00564
00565 pHit0 = aMucRecHitContainer->GetHit(iPart, iSeg, iGap0, iHit0);
00566 if (!pHit0) {
00567 log << MSG::FATAL << "MucRecRoadFinder-E10 " << " part " << iPart
00568 << " seg " << iSeg << " gap " << iGap0 << " null pointer"
00569 << endl;
00570 }
00571 else {
00572
00573
00574 if(m_united == 1 && pHit0->GetHitMode() != -1) continue;
00575
00576 count1 = aMucRecHitContainer->GetGapHitCount(iPart, iSeg, iGap1);
00577 if(m_MsOutput) cout << "part " << iPart << " seg " << iSeg << " orient " << iOrient
00578 << " gap0 " << iGap0 << " count0 " << count0
00579 << " gap1 " << iGap1 << " count1 " << count1 << endl;
00580 for (int iHit1 = 0; iHit1 < count1; iHit1++) {
00581
00582 pHit1 = aMucRecHitContainer->GetHit(iPart, iSeg, iGap1, iHit1);
00583 if (!pHit1) {
00584 log << MSG::ERROR << "MucRecRoadFinder-E10 " << " part " << iPart
00585 << " seg " << iSeg << " gap " << iGap1 << " null pointer"
00586 << endl;
00587 } else {
00588
00589
00590 if(m_united == 1 && pHit1->GetHitMode() != -1) continue;
00591
00592
00593
00594 MucRec2DRoad *road2D = new MucRec2DRoad(iPart, iSeg, iOrient, 0.0, 0.0, 3000.0);
00595 if (!road2D) {
00596 log << MSG::FATAL << "MucRecRoadFinder-E20 " << "failed to create 2D road!" << endreq;
00597 continue;
00598 }
00599 road2D->SetMaxNSkippedGaps(kMaxSkippedGaps);
00600
00601 if (!pHit0->GetGap()) log << MSG::ERROR << "pHit0->GetGap(), null pointer" << endreq;
00602 if (pHit0->GetGap()->Orient() == iOrient) {
00603
00604
00605
00606
00607
00608 bool isseed = true;
00609 pHit0->SetHitSeed(true);
00610 pHit1->SetHitSeed(true);
00611
00612 road2D->AttachHit(pHit0);
00613 if(m_MsOutput) cout << "pHit0 attached, " << road2D->GetTotalHits()
00614 << ", hitId= "<<pHit0->Part()<<" "<<pHit0->Seg()<<" "<<pHit0->Gap()<<" "<<pHit0->Strip()<<endl;
00615 }
00616
00617 if (pHit1->GetGap()->Orient() == iOrient) {
00618
00619
00620
00621
00622 road2D->AttachHit(pHit1);
00623 if(m_MsOutput) cout << "pHit1 attached, " << road2D->GetTotalHits()
00624 << ", hitId= "<<pHit1->Part()<<" "<<pHit1->Seg()<<" "<<pHit1->Gap()<<" "<<pHit1->Strip()<<endl;
00625 }
00626 if(m_MsOutput) cout << "Hit cout " << road2D->GetTotalHits() << ", slope " << road2D->GetSlope() << endl;
00627
00628
00629
00630
00631
00632
00633 for (int i = 0; i < length; i++) {
00634 int iGap = kSearchOrder[iLoop][iPart][iOrient][i];
00635
00636 float dXmin = kInfinity;
00637 MucRecHit *pHit = 0;
00638
00639 float Win = kWindowSize[iPart][iGap];
00640
00641
00642
00643
00644
00645
00646
00647 count = aMucRecHitContainer->GetGapHitCount(iPart, iSeg, iGap);
00648 for (int iHit = 0; iHit < count; iHit++) {
00649 pHit = aMucRecHitContainer->GetHit(iPart, iSeg, iGap, iHit);
00650 if (!pHit) {
00651 log << MSG::FATAL << "MucRecRoadFinder-E20 null cluster pointer pHit" << endreq;
00652 continue;
00653 }
00654 else {
00655
00656
00657 if(m_united == 1 && pHit->GetHitMode() != -1) continue;
00658
00659
00660 float dX = road2D->GetHitDistance(pHit);
00661 if(m_MsOutput) cout<<"dX = "<<dX<<" Win = "<<Win<<", hit: "<<pHit->Part()<<" "<<pHit->Seg()<<" "<<pHit->Gap()<<" "<<pHit->Strip()<<endl;
00662
00663 if (dX < Win) {
00664
00665
00666 if(iGap == iGap0 || iGap == iGap1) {
00667 if(pHit->GetHitMode() == -1){
00668 pHit->SetHitMode(3);
00669 }
00670 else if(pHit->GetHitMode() == 3){
00671 pHit->SetHitMode(4);
00672 }
00673 }
00674
00675 if(m_onlyseedfit == 0)road2D->AttachHit(pHit);
00676 else {
00677 if(iGap == iGap0 || iGap == iGap1) road2D->AttachHit(pHit);
00678 else road2D->AttachHitNoFit(pHit);
00679 }
00680 }
00681 }
00682 }
00683
00684
00685 }
00686
00687
00688
00689
00690
00691
00692
00693
00694 bool lastGapOK = false;
00695 if ( road2D->GetLastGap() >= kMinLastGap2D ) {
00696 lastGapOK = true;
00697 }
00698 else {
00699 log<<MSG::INFO << " 2D lastGap " << road2D->GetLastGap() << " < " << kMinLastGap2D << endreq;
00700 }
00701
00702
00703
00704
00705 bool firedGapsOK = false;
00706 if (road2D->GetNGapsWithHits() >= kMinFiredGaps) {
00707 firedGapsOK = true;
00708 }
00709 else {
00710 log<<MSG::INFO << " 2D firedGaps " << road2D->GetNGapsWithHits() << " < " << kMinFiredGaps << endreq;
00711 }
00712
00713
00714 bool duplicateRoad = false;
00715 int maxSharedHits = 0, sharedHits = 0;
00716 if (iOrient == 0) {
00717 for (int index = 0; index < (int)p2DRoad0C->size(); index++) {
00718 sharedHits =(*p2DRoad0C)[index]->GetNSharedHits(road2D);
00719 }
00720 if (maxSharedHits < sharedHits) maxSharedHits = sharedHits;
00721 }
00722 else {
00723 for (int index = 0; index < (int)p2DRoad1C->size(); index++) {
00724 sharedHits = (*p2DRoad1C)[index]->GetNSharedHits(road2D);
00725 }
00726 if (maxSharedHits < sharedHits) maxSharedHits = sharedHits;
00727 }
00728
00729 if(road2D->GetTotalHits() == maxSharedHits
00730 || maxSharedHits >= kMinSharedHits2D ) {
00731 duplicateRoad = true;
00732 log<<MSG::INFO << " maxSharedHits " << maxSharedHits << " > " << kMinSharedHits2D << endreq;
00733 }
00734
00735
00736
00737
00738
00739
00740
00741 if (lastGapOK && firedGapsOK && !duplicateRoad) {
00742 if (iOrient == 0) {
00743 log<<MSG::INFO << "Add a road0" << endreq;
00744 p2DRoad0C->add(road2D);
00745 }
00746 else {
00747 log<<MSG::INFO << "Add a road1" << endreq;
00748 p2DRoad1C->add(road2D);
00749 }
00750 }
00751 else {
00752
00753
00754 vector<MucRecHit*> road2DHits = road2D->GetHits();
00755 for(int i=0; i< road2DHits.size(); i++)
00756 {
00757 MucRecHit *ihit = road2DHits[i];
00758 if(ihit->Gap() == iGap0 || ihit->Gap() == iGap1){
00759 if(ihit->GetHitMode() == 3) ihit->SetHitMode(-1);
00760 if(ihit->GetHitMode() == 4) ihit->SetHitMode(3);
00761 }
00762 }
00763 delete road2D;
00764 }
00765 }
00766 }
00767 }
00768 }
00769 }
00770 }
00771 }
00772 }
00773
00774 int print2DRoadInfo = 0;
00775 if (print2DRoadInfo == 1) {
00776
00777 cout << "In 2DRoad container : " << endl
00778 << " Num of 2DRoad0 = " << p2DRoad0C->size() << endl
00779 << endl;
00780
00781 for (int iRoad = 0; iRoad < (int)p2DRoad0C->size(); iRoad++) {
00782
00783 MucRec2DRoad *road = (*p2DRoad0C)[iRoad];
00784 cout << " " << iRoad << "th 2DRoad0 :" << endl
00785 << " Part = " << road->GetPart() << endl
00786 << " Seg = " << road->GetSeg() << endl
00787 << " Orient = " << road->GetOrient() << endl
00788 << " LastGap = " << road->GetLastGap() << endl
00789 << " TotalHits = " << road->GetTotalHits() << endl
00790 << " DOF = " << road->GetDegreesOfFreedom() << endl
00791 << " Slope = " << road->GetSlope() << endl
00792 << " Intercept = " << road->GetIntercept() << endl
00793 << endl;
00794
00795 }
00796
00797 cout << " Num of 2DRoad1 = " << p2DRoad1C->size() << endl
00798 << endl;
00799
00800 for ( int iRoad = 0; iRoad < (int)p2DRoad1C->size(); iRoad++) {
00801
00802 MucRec2DRoad *road = (*p2DRoad1C)[iRoad];
00803 cout << " " << iRoad << "th 2DRoad1 :" << endl
00804 << " Part = " << road->GetPart() << endl
00805 << " Seg = " << road->GetSeg() << endl
00806 << " Orient = " << road->GetOrient() << endl
00807 << " LastGap = " << road->GetLastGap() << endl
00808 << " TotalHits = " << road->GetTotalHits() << endl
00809 << " DOF = " << road->GetDegreesOfFreedom() << endl
00810 << " Slope = " << road->GetSlope() << endl
00811 << " Intercept = " << road->GetIntercept() << endl
00812 << endl;
00813
00814 }
00815 }
00816
00817
00818
00819 for ( int iRoad0 = 0; iRoad0 < (int)p2DRoad0C->size(); iRoad0++) {
00820 MucRec2DRoad *road0 = (*p2DRoad0C)[iRoad0];
00821 for ( int iRoad1 = 0; iRoad1 < (int)p2DRoad1C->size(); iRoad1++){
00822 MucRec2DRoad *road1 = (*p2DRoad1C)[iRoad1];
00823
00824
00825 if ( !(road0 &&road1) ) {
00826 cout << "Null pointer to road0 or road1: "
00827 << "road0 = " << road0
00828 << "road1 = " << road1
00829 << endl;
00830 continue;
00831 }
00832
00833
00834 if ( (road0->GetPart() != road1->GetPart())
00835 || (road0->GetSeg() != road1->GetSeg()) ) {
00836 continue;
00837 }
00838
00839 MucRec3DRoad *road3D = new MucRec3DRoad(road0, road1);
00840
00841
00842
00843
00844
00845
00846
00847
00848
00849
00850
00851
00852 bool lastGapDeltaOK = false;
00853 if ( road3D->GetLastGapDelta() <= kMaxDeltaLastGap ) {
00854 lastGapDeltaOK = true;
00855 }
00856 else {
00857 cout << "LastGapDelta = " << road3D->GetLastGapDelta() << endl;
00858 }
00859
00860 bool totalHitsDeltaOK = false;
00861 if ( road3D->GetTotalHitsDelta() <= kMaxDeltaTotalHits ) {
00862 totalHitsDeltaOK = true;
00863 }
00864 else {
00865
00866 }
00867
00868 bool chiSquareCutOK = false;
00869 float xChiSquare = road0->GetReducedChiSquare();
00870 float yChiSquare = road1->GetReducedChiSquare();
00871 if ( xChiSquare <= kMaxXChisq
00872 && yChiSquare <= kMaxYChisq ) {
00873 chiSquareCutOK = true;
00874 }
00875 else {
00876 cout << "xChiSquare = " << xChiSquare
00877 << "yChiSquare = " << yChiSquare
00878 << endl;
00879 }
00880
00881 bool minLastGapOK = false;
00882 if ( road3D->GetLastGap() >= kMinLastGap3D ) {
00883 minLastGapOK = true;
00884 }
00885 else {
00886 log<<MSG::INFO << " minLastGap = " << road3D->GetLastGap() << endreq;
00887 }
00888
00889 bool duplicateRoad = false;
00890 int maxSharedHits = 0, sharedHits = 0;
00891 for ( int i = 0; i < (int)p3DRoadC->size(); i++){
00892 sharedHits =(*p3DRoadC)[i]->GetNSharedHits(road3D);
00893 if ( maxSharedHits < sharedHits) maxSharedHits = sharedHits;
00894
00895
00896 }
00897 if(road3D->GetTotalHits() == maxSharedHits
00898 || maxSharedHits >= kMinSharedHits2D ) {
00899 duplicateRoad = true;
00900 log<<MSG::INFO << " MaxShareHits = " << maxSharedHits << endreq;
00901 }
00902
00903 if ( lastGapDeltaOK &&
00904 totalHitsDeltaOK &&
00905 chiSquareCutOK &&
00906 minLastGapOK &&
00907 !duplicateRoad ) {
00908 float vx, vy, x0, y0;
00909 float slope1 = 100, slope0 = 100;
00910 if(fabs(road1->GetSlope())<100) slope1 = road1->GetSlope();
00911 if(fabs(road0->GetSlope())<100) slope0 = road0->GetSlope();
00912
00913 if ( road3D->GetPart() == 1 ){
00914 road3D->TransformPhiRToXY( slope1, slope0,
00915 road1->GetIntercept(), road0->GetIntercept(),
00916 vx, vy, x0, y0);
00917 }
00918 else {
00919 vx = slope1;
00920 x0 = road1->GetIntercept();
00921 vy = slope0;
00922 y0 = road0->GetIntercept();
00923 }
00924 road3D->SetSimpleFitParams(vx, x0, vy, y0);
00925
00926 log<<MSG::INFO << "Add a 3D Road ... " << endreq;
00927
00928 float startx = 0.0, starty = 0.0, startz = 0.0;
00929 float sigmax = 0.0, sigmay = 0.0, sigmaz = 0.0;
00930 road3D->ProjectWithSigma(0, startx, starty, startz, sigmax, sigmay, sigmaz);
00931
00932
00933
00934
00935 float vz = 1;
00936 float sign = vy/fabs(vy);
00937 float signx = vx/fabs(vx);
00938
00939 if(road3D->GetPart() == 1){
00940 if(road3D->GetSeg() > 4){
00941
00942 vx *= -sign;
00943 vy *= -sign;
00944 vz *= -sign;
00945 }
00946 else if(road3D->GetSeg()<2){
00947 vx *= signx;
00948 vy *= signx;
00949 vz *= signx;
00950 }
00951 else if(road3D->GetSeg()>=3 && road3D->GetSeg()<=4){
00952 vx *= -signx;
00953 vy *= -signx;
00954 vz *= -signx;
00955 }
00956 else{
00957 vx *= sign;
00958 vy *= sign;
00959 vz *= sign;
00960 }
00961 }
00962 else if(road3D->GetPart() == 0){
00963
00964
00965
00966
00967
00968 }
00969 else if(road3D->GetPart() == 2){
00970
00971
00972 vx *= -1;
00973 vy *= -1;
00974 vz *= -1;
00975
00976
00977
00978
00979 }
00980
00981
00982 Hep3Vector mom(vx, vy, vz);
00983
00985
00986
00987
00988
00989 startx /= 10; starty /= 10; startz /= 10;
00990 startx -= vx/mom.mag(); starty -= vy/mom.mag(); startz -= vz/mom.mag();
00991
00992
00993 RecMucTrack *aTrack = new RecMucTrack();
00994 aTrack->SetExtMucPos(startx, starty, startz);
00995 aTrack->SetExtMucMomentum(vx, vy, vz);
00996 aTrack->SetMucPos(startx, starty, startz);
00997 aTrack->SetMucMomentum(vx, vy, vz);
00998 aTrack->SetCurrentPos( startx, starty, startz);
00999 aTrack->SetCurrentDir( vx, vy, vz);
01000 aTrack->SetRecMode(3);
01001
01002
01003
01004 aTrack->setTrackId(trackId);
01005 aTrack->setId(muctrackId);
01006 trackId++;
01007 muctrackId++;
01008
01009
01010 aMucTrackCol->add(aTrack);
01011 TrackFinding(aTrack);
01012 p3DRoadC->add(road3D);
01013
01015
01016 vector<MucRecHit*> attachedHits = aTrack->GetHits();
01017 vector<MucRecHit*> expectedHits = aTrack->GetExpectedHits();
01018 vector<float> distanceHits = aTrack->getDistHits();
01019
01020 for(int i=0; i< expectedHits.size(); i++)
01021 {
01022 MucRecHit *ihit = expectedHits[i];
01023
01024 }
01025
01026 vector<int> multiHit;
01027 for(int i=0; i< attachedHits.size(); i++)
01028 {
01029 MucRecHit *ihit = attachedHits[i];
01030
01031
01032 int hitnum = 0;
01033 for(int k=0; k < attachedHits.size(); k++){
01034 MucRecHit *khit = attachedHits[k];
01035 if((ihit->Part()==khit->Part())&&(ihit->Seg()==khit->Seg())&&(ihit->Gap()==khit->Gap()))
01036 hitnum++;
01037 }
01038 multiHit.push_back(hitnum);
01039
01040
01041 }
01042
01043 for(int i=0; i< expectedHits.size(); i++)
01044 {
01045
01046 MucRecHit *ihit = expectedHits[i];
01047
01048
01049 int diff = -999;
01050
01051 for(int j=0; j< attachedHits.size(); j++){
01052 MucRecHit *jhit = attachedHits[j];
01053
01054
01055
01056 if((jhit->Part()==ihit->Part())&&(jhit->Seg()==ihit->Seg())&&(jhit->Gap()==ihit->Gap())&&attachedHits.size()==distanceHits.size())
01057 {
01058 diff = ihit->Strip() - jhit->Strip();
01059
01060
01061 if(m_NtOutput>=2){
01062
01063 m_part = ihit->Part(); m_seg = ihit->Seg(); m_gap = ihit->Gap();
01064 m_strip = jhit->Strip();
01065 m_diff = diff;
01066 m_dist = distanceHits[j];
01067
01068
01069 m_angle_cosmic = -999;
01070 m_angle_updown = -999;
01071
01072
01073
01074 m_ngapwithhits = aTrack->numLayers();
01075 m_nhit = aTrack->numHits();
01076 m_maxhit = aTrack->maxHitsInLayer();
01077 m_multihit = multiHit[j];
01078 m_run = eventHeader->runNumber();
01079 m_event = eventHeader->eventNumber();
01080
01081 m_tuple->write();
01082 }
01083 }
01084
01085
01086 }
01087
01088
01089
01090 if(diff == -999) {
01091 if(m_NtOutput>=2){
01092 m_part = ihit->Part(); m_seg = ihit->Seg(); m_gap = ihit->Gap();
01093 m_strip = ihit->Strip();
01094 m_diff = diff;
01095 m_dist = -999;
01096 m_angle_updown = -999;
01097 m_angle_cosmic = -999;
01098
01099
01100
01101 m_ngapwithhits = aTrack->numLayers();
01102 m_run = eventHeader->runNumber();
01103 m_event = eventHeader->eventNumber();
01104
01105 }
01106 }
01107
01108
01109 }
01111
01112
01113
01114
01115
01116
01117
01118
01119
01120
01121
01122
01123
01124
01125
01126
01127
01128
01129
01130
01131
01132
01133
01134
01135
01136
01137
01138
01139
01141
01142
01143
01144 }
01145 else {
01146
01147 delete road3D;
01148
01149 }
01150 }
01151 }
01152
01153
01154
01155 RecMucTrack *aaTrack = 0;
01156 RecMucTrack *bbTrack = 0;
01157
01158 int hasMucUp = 0;
01159 int hasMucDown = 0;
01160 for(int iTrack = 0; iTrack < (int)aMucTrackCol->size(); iTrack++) {
01161 aaTrack = (*aMucTrackCol)[iTrack];
01162 if((aaTrack->getMucMomentum()).phi()<3.1415926&&(aaTrack->getMucMomentum()).phi()>0) hasMucUp++;
01163 else hasMucDown++;
01164
01165
01166 double px,py,pz,p,theta,phi;
01167 px = (aaTrack->getMucMomentum()).x();
01168 py = (aaTrack->getMucMomentum()).y();
01169 pz = (aaTrack->getMucMomentum()).z();
01170
01171 if(py<0)continue;
01172
01173 if(m_NtOutput>=1){
01174
01175 m_angle_updown = -999;
01176 m_angle_cosmic = cosmicMom.angle(aaTrack->getMucMomentum());
01177 if(m_angle_cosmic>1.57) m_angle_cosmic = 3.14159265 - m_angle_cosmic;
01178 m_run = eventHeader->runNumber();
01179 m_event = eventHeader->eventNumber();
01180
01181 Hep3Vector rotate(-px,pz,py);
01182 theta = rotate.theta();
01183 phi = rotate.phi();
01184
01185 m_px = px; m_py = py; m_pz = pz;
01186 m_theta = theta; m_phi = phi;
01187 m_theta_pipe = (aaTrack->getMucMomentum()).theta();
01188 m_phi_pipe = (aaTrack->getMucMomentum()).phi();
01189
01190
01191
01192
01193 Hep3Vector mucPos = aaTrack->getMucPos();
01194 double posx, posy, posz;
01195 posx = mucPos.x(); posy = mucPos.y(); posz = mucPos.z();
01196
01197 m_projx = -999; m_projz = -999;
01198 if(py!=0){
01199 m_projx = posx - px/py*posy;
01200 m_projz = posz - pz/py*posy;
01201 }
01202
01203 }
01204
01205 }
01206 if(m_NtOutput>=1){
01207 m_mucUp = hasMucUp; m_mucDown = hasMucDown;
01208 m_tuple->write();
01209 }
01210
01211 int forCosmic = 0;
01212 if(aMucTrackCol->size()>=2 && forCosmic == 1){
01213 for(int iTrack = 0; iTrack < (int)aMucTrackCol->size(); iTrack++) {
01214 log << MSG::DEBUG << "iTrack " << iTrack << " / " <<(int)aMucTrackCol->size()<<endreq;
01215 aaTrack = (*aMucTrackCol)[iTrack];
01216 if(aaTrack->trackId()>=0) continue;
01217 aaTrack->setTrackId(iTrack);
01218
01219 for(int jTrack = iTrack+1; jTrack < (int)aMucTrackCol->size(); jTrack++){
01220 bbTrack = (*aMucTrackCol)[jTrack];
01221
01222 Hep3Vector mom_a = aaTrack->getMucMomentum();
01223 Hep3Vector mom_b = bbTrack->getMucMomentum();
01224
01225
01226 if(fabs(mom_a.angle(mom_b) - 3.1415926)<0.2)
01227 {
01228 bbTrack->setTrackId(iTrack);
01229
01230
01231
01232
01233
01234 }
01235
01236 if(m_NtOutput>=1){
01237 m_angle_cosmic = cosmicMom.angle(mom_a);
01238 if(cosmicMom.angle(mom_a)>1.57) m_angle_cosmic = 3.14159265 - cosmicMom.angle(mom_a);
01239
01240 m_angle_updown = fabs(mom_a.angle(mom_b) - 3.1415926);
01241 m_px = -999; m_py = -999; m_pz = -999; m_theta = -999; m_phi = -999; m_theta_pipe = -999; m_phi_pipe = -999;
01242 m_px_mc = -999; m_py_mc = -999; m_pz_mc = -999; m_theta_mc = -999; m_phi_mc = -999;
01243 m_theta_mc_pipe = -999; m_phi_mc_pipe = -999;
01244
01245
01246 }
01247 }
01248
01249 }
01250
01251
01252 }
01253
01254 if( p3DRoadC->size() !=0 ) {
01255 log<<MSG::INFO << "In 3DRoad container : "
01256 << " Num of 3DRoad = " << p3DRoadC->size() <<endreq;
01257
01258 int print2DRoadInfo = 0;
01259 if (print2DRoadInfo == 1) {
01260 for ( int iRoad = 0; iRoad < (int)p3DRoadC->size(); iRoad++) {
01261 MucRec3DRoad *road = (*p3DRoadC)[iRoad];
01262 cout << endl << " " << iRoad << "th 3DRoad :" << endl
01263 << " Part = " << road->GetPart() << endl
01264 << " Seg = " << road->GetSeg() << endl
01265 << " NumGapsWithHits = " << road->GetNGapsWithHits() << endl
01266 << " TotalHits = " << road->GetTotalHits() << endl
01267 << " MaxHitsPerGap = " << road->GetMaxHitsPerGap() << endl
01268 << " LastGapDelta = " << road->GetLastGapDelta() << endl
01269 << " TotalHitsDelta = " << road->GetTotalHitsDelta() << endl
01270 << " DegreesOfFreedom = " << road->GetDegreesOfFreedom() << endl
01271 << " ReducedChiSquare = " << road->GetReducedChiSquare() << endl
01272 << " SlopeZX = " << road->GetSlopeZX() << endl
01273 << " InterceptZX = " << road->GetInterceptZX() << endl
01274 << " SlopeZY = " << road->GetSlopeZY() << endl
01275 << " InterceptZY = " << road->GetInterceptZY() << endl
01276 << " Hits Info : " << endl;
01277
01278 }
01279 }
01280
01281 m_NEventReco ++;
01282 }
01283
01284
01285 for(int i = 0 ; i < p3DRoadC->size(); i++)
01286 delete (*p3DRoadC)[i];
01287 for(int i = 0 ; i < p2DRoad0C->size(); i++)
01288 delete (*p2DRoad0C)[i];
01289 for(int i = 0 ; i < p2DRoad1C->size(); i++)
01290 delete (*p2DRoad1C)[i];
01291
01292 p3DRoadC->clear();
01293 p2DRoad0C->clear();
01294 p2DRoad1C->clear();
01295 delete p3DRoadC;
01296 delete p2DRoad0C;
01297 delete p2DRoad1C;
01298 return StatusCode::SUCCESS;
01299 }