MucRecRoadFinder Class Reference

#include <MucRecRoadFinder.h>

List of all members.

Public Member Functions

 MucRecRoadFinder (const std::string &name, ISvcLocator *pSvcLocator)
 ~MucRecRoadFinder ()
StatusCode initialize ()
StatusCode execute ()
StatusCode finalize ()
void TrackFinding (RecMucTrack *aTrack)

Private Attributes

int m_NEvent
int m_NEventWithHit
int m_NEventReco
int m_NHitsLostTotal
int m_NHitsTotal
std::vector< int > m_NHitsLost
std::vector< int > m_NHitsLostInGap
int m_fittingMethod
std::string m_configFile
int m_mccosmic
int m_NtOutput
int m_onlyseedfit
int m_maxHitsRec
int m_united
int m_seedtype
bool m_MsOutput
std::string m_filter_filename
std::vector< FilterEventm_filter_event
MucRecHitContaineraMucRecHitContainer
NTuple::Tuple * m_tuple
NTuple::Item< double > m_part
NTuple::Item< double > m_seg
NTuple::Item< double > m_gap
NTuple::Item< double > m_strip
NTuple::Item< double > m_diff
NTuple::Item< double > m_dist
NTuple::Item< double > m_run
NTuple::Item< double > m_event
NTuple::Item< double > m_ngapwithhits
NTuple::Item< double > m_nhit
NTuple::Item< double > m_maxhit
NTuple::Item< double > m_multihit
NTuple::Item< double > m_angle_cosmic
NTuple::Item< double > m_angle_updown
NTuple::Item< double > m_px
NTuple::Item< double > m_py
NTuple::Item< double > m_pz
NTuple::Item< double > m_theta
NTuple::Item< double > m_phi
NTuple::Item< double > m_theta_pipe
NTuple::Item< double > m_phi_pipe
NTuple::Item< double > m_px_mc
NTuple::Item< double > m_py_mc
NTuple::Item< double > m_pz_mc
NTuple::Item< double > m_theta_mc
NTuple::Item< double > m_phi_mc
NTuple::Item< double > m_theta_mc_pipe
NTuple::Item< double > m_phi_mc_pipe
NTuple::Item< double > m_emcUp
NTuple::Item< double > m_emcDown
NTuple::Item< double > m_mucUp
NTuple::Item< double > m_mucDown
NTuple::Item< double > m_projx
NTuple::Item< double > m_projz

Classes

struct  FilterEvent


Detailed Description

Reconstruction of Muon Chamber tracks by combining two 2DRoads to form a 3DRoad.

Author:
Zhengyun You {mailto:youzy@hep.pku.cn}

Definition at line 28 of file MucRecRoadFinder.h.


Constructor & Destructor Documentation

MucRecRoadFinder::MucRecRoadFinder ( const std::string name,
ISvcLocator *  pSvcLocator 
)

Definition at line 49 of file MucRecRoadFinder.cxx.

References m_configFile, m_filter_filename, m_fittingMethod, m_maxHitsRec, m_mccosmic, m_MsOutput, m_NtOutput, m_onlyseedfit, m_seedtype, and m_united.

00049                                                                                   :
00050   Algorithm(name, pSvcLocator)
00051 {
00052   // Declare the properties  
00053   declareProperty("FittingMethod", m_fittingMethod = 2);
00054   declareProperty("ConfigFile",  m_configFile = "MucConfig.xml");
00055   declareProperty("McCosmic",    m_mccosmic = 0);
00056   declareProperty("OnlySeedFit",    m_onlyseedfit = 0);
00057   declareProperty("NtOutput",      m_NtOutput = 0);
00058   declareProperty("MaxHitsRec",      m_maxHitsRec = 200);
00059   declareProperty("United", m_united = 0);
00060   declareProperty("SeedType", m_seedtype = 0);
00061   declareProperty("MsOutput", m_MsOutput = false);
00062   declareProperty("FilterFile", m_filter_filename);
00063 }

MucRecRoadFinder::~MucRecRoadFinder (  )  [inline]

Definition at line 32 of file MucRecRoadFinder.h.

00032 {};


Member Function Documentation

StatusCode MucRecRoadFinder::execute (  ) 

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   // Part 1: Get the event header, print out event and run number
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 //cout<<"events  run  "<<event <<"  "<<run<<endl;
00190   string release = getenv("BES_RELEASE");
00191   // if(release=="6.6.5"&&run==35896&&event==7448){return ( StatusCode::SUCCESS);}
00192   // if(release=="6.6.5"&&run==37241&&event==1690731){return ( StatusCode::SUCCESS);}
00193   // filter the event
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   //Part 2: Retrieve MC truth 
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       //cout<<"in mucroadfinder pid = "<<pid<<endl;
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         //m_tuple->write();  
00242       }
00243 
00244       //if(initialMomentum.py()>0)cout<<"py>0?????  "<<pid<<endl;
00245       cosmicMom = startP;
00246 
00247     }
00248   }
00249   /*
00250   //Part 3: Retrieve MDC digi 
00251   SmartDataPtr<MdcDigiCol> mdcDigiCol(eventSvc(),"/Event/Digi/MdcDigiCol");
00252   if (!mdcDigiCol) {
00253   log << MSG::FATAL << "Could not find MDC digi" << endreq;
00254   return( StatusCode::FAILURE);
00255   }
00256 
00257   MdcDigiCol::iterator iter1 = mdcDigiCol->begin();
00258   digiId = 0; 
00259 
00260   for (;iter1 != mdcDigiCol->end(); iter1++, digiId++) {
00261   log << MSG::INFO << "MDC digit No: " << digiId << endreq;
00262 
00263   log << MSG::INFO 
00264   << " time_channel = " << (*iter1)->getTimeChannel()
00265   << " charge_channel = " << (*iter1)->getChargeChannel() 
00266   << endreq;   
00267   }
00268 
00269   //Part 4: Retrieve TOF digi 
00270   SmartDataPtr<TofDigiCol> tofDigiCol(eventSvc(),"/Event/Digi/TofDigiCol");
00271   if (!tofDigiCol) {
00272   log << MSG::FATAL << "Could not find TOF digi" << endreq;
00273   return( StatusCode::FAILURE);
00274   }
00275 
00276   TofDigiCol::iterator iter2 = tofDigiCol->begin();
00277   digiId = 0; 
00278 
00279   for (;iter2 != tofDigiCol->end(); iter2++, digiId++) {
00280   log << MSG::INFO << "TOF digit No: " << digiId << endreq;
00281 
00282   }
00283 
00284   //Part 5: Retrieve EMC digi 
00285   SmartDataPtr<EmcDigiCol> emcDigiCol(eventSvc(),"/Event/Digi/EmcDigiCol");
00286   if (!emcDigiCol) {
00287   log << MSG::FATAL << "Could not find EMC digi" << endreq;
00288   return( StatusCode::FAILURE);
00289   }
00290 
00291   EmcDigiCol::iterator iter3 = emcDigiCol->begin();
00292   digiId = 0; 
00293 
00294   for (;iter3 != emcDigiCol->end(); iter3++, digiId++) {
00295   log << MSG::INFO << "Emc digit No: " << digiId << endreq;
00296 
00297   log << MSG::INFO 
00298   << " time_channel = " << (*iter3)->getTimeChannel()
00299   << " charge_channel = " << (*iter3)->getChargeChannel() 
00300   << endreq;   
00301   }
00302    */
00303 
00304   //Part 6: Retrieve MUC digi 
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)  // if this algorithm run individualy, we need create hitcollection and trackcollection...
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       //aMucRecHitContainer->AddHit(mucID);
00354       //if(part==1 && seg==2 && gap==8 && (strip==19||strip==16));//cout<<"noise hit!!!=========="<<endl;
00355       //else            /////this problem has been solved!
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       //then register Recon
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     // Create track Container;
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     // check MucTrackCol registered
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;  //too many digit! abnormal------2008-04-17 
00410     //********************************//
00411     // Create track Container;
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     // Retrieve Ext track Col from TDS for setting trackId 
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     //return( StatusCode::FAILURE);
00460   }
00461   else{
00462     RecEmcShowerCol::iterator iShowerCol;
00463     for(iShowerCol=aShowerCol->begin();
00464         iShowerCol!=aShowerCol->end();
00465         iShowerCol++){
00466       /*
00467          cout<<"check EmcRecShowerCol registered:"<<endl;
00468          cout<<"shower id: "<<(*iShowerCol)->ShowerId()<<"\t"
00469          <<"shower energy: "<<(*iShowerCol)->Energy()<<"\n"
00470          <<"shower position: "<<(*iShowerCol)->Position()<<endl;
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   // Find or create the 3D road container ...
00484   if (!p3DRoadC) {
00485     p3DRoadC = new MucRec3DRoadContainer();
00486   }
00487   p3DRoadC->clear();    // make sure that the container is empty
00488 
00489   MucRec2DRoadContainer *p2DRoad0C = new MucRec2DRoadContainer();
00490   // Find or create the 2D road0 container ...
00491   if (!p2DRoad0C) {
00492     p2DRoad0C = new MucRec2DRoadContainer();
00493   }
00494   p2DRoad0C->clear();    // make sure that the container is empty
00495 
00496   MucRec2DRoadContainer *p2DRoad1C = new MucRec2DRoadContainer();
00497   // Find or create the 2D road1 container ...
00498   if (!p2DRoad1C) {
00499     p2DRoad1C = new MucRec2DRoadContainer();
00500   }
00501   p2DRoad1C->clear();    // make sure that the container is empty
00502   log << MSG::INFO << "Muc 2D 3D Road Container created" << endreq;
00503 
00504   // We have all of the input and output parameters now, so do
00505   // something useful with them!
00506   //static bool  first_call = true;
00507   //
00508   // Now find some roads!
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++) { //now only 1 loop
00521           // checkout the seed gap(s) from search order
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 {  //new method to calc seed gaps------//
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               }//make sure that iGap0 record the 1st gap in this orient
00546               if(part == iPart && seg == iSeg && orient == iOrient && setgap0 == 1 && gap != iGap0 ){
00547                 iGap1 = gap;
00548               }//make sure that iGap1 record the last gap in this orient
00549             }     
00550           }
00551           //                                                              
00552           if(m_MsOutput) cout <<"Find seed gap in ( "<<iPart<<"  "<<iSeg<<" ) gap0: "<<iGap0<<"  gap1: "<<iGap1<<endl;
00553 
00554           //new method to calc seed gaps------//
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             //cout << "iHit0 = " << iHit0 << endl;
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               //this algo use with ext and this hit has been used before; 
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                 //cout << "iHit1 = " << iHit1 << endl;
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                   //this algo use with ext and this hit has been used before;                   
00590                   if(m_united == 1 && pHit1->GetHitMode() != -1) continue;
00591 
00592                   // Found seed hit(s), create a new 2D road object,
00593                   // and attach hit pHit1 and pHit0
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                     //float xStrip, yStrip, zStrip;
00604                     //MucRecGeometry::GetPointer()->GetStripPointer(pHit0->GetSoftID())->GetCenterPos(xStrip, yStrip, zStrip);
00605                     //cout << " x = " << xStrip << " y = " << yStrip << " z = " << zStrip << endl;
00606 
00607                     //set hit as seed then
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                     //cout << pHit1->GetSoftID() << endl;
00619                     //float xStrip, yStrip, zStrip;
00620                     //MucRecGeometry::GetPointer()->GetStripPointer(pHit1->GetSoftID())->GetCenterPos(xStrip, yStrip, zStrip);
00621                     //cout << " x = " << xStrip << " y = " << yStrip << " z = " << zStrip << endl;
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                   // After first (two) hit(s) is(are) attached,
00629                   // calculate reference pos, direction.
00630                   // Project to the other planes;
00631                   // attach clusters within search window.
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                     //Win = road2D->GetSearchWindowSize(iGap);  
00641 
00642                     // following line should be removed after
00643                     // we have a good function classes for search window
00644                     //Win = WindowParGamma[iGap];
00645 
00646                     //search hit in iGap
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                         //this algo use with ext and this hit has been used before;                   
00657                         if(m_united == 1 && pHit->GetHitMode() != -1) continue;
00658 
00659                         // Get the displacement of hit pHit to road2D
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                           // Attach the hit if it exists
00665                           // cout << "meet window " << pHit->GetSoftID() << endl;
00666                           if(iGap == iGap0 || iGap == iGap1) { //hits in these gap is used to be seed, so unvailable for calibrate effi;
00667                             if(pHit->GetHitMode() == -1){  
00668                               pHit->SetHitMode(3); //self road finder mode
00669                             }
00670                             else if(pHit->GetHitMode() == 3){
00671                               pHit->SetHitMode(4); //self road finder mode used!!!
00672                             }
00673                           }
00674                           // attach hits if we need fit or hit. OR, only attach hits not in seed layers
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                   } //  for (int i = 0; i < length; ...), Search all gaps
00686 
00687 
00688                   // Ok we found a road already, we need to know
00689                   // whether we should save it or not
00690                   // check road quality...
00691                   //
00692 
00693                   // 1. last gap of the road
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                   // 2. positon at reference plane
00703 
00704                   // 3. number of gaps contain hits
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                   // 4. duplicated road check
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                   // Here is our criteria for a candidate road
00736                   // i.e. 1) There are at least two gaps containing hits
00737                   //      2) LastGap of the road passes last plane cut
00738                   //      3) Reference position passes vertex cut
00739                   //      4) No existing road share same hits with the road
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                     // reset hit mode to be -1 if this road useless! 
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                 } // else { 
00766               } // for ( int iHit = 0 ...) 
00767               } // else {
00768             } // for ( int iHit0 ...)
00769             } // for (int iLoop ...)
00770           } // for (int iOrient ...)
00771         } // for (int iSeg ...)
00772       } // for(iPartArm ...)
00773 
00774       int print2DRoadInfo = 0;
00775       if (print2DRoadInfo == 1) {
00776         // print 2DRoad container info.
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           //road->PrintHitsInfo();
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           //road->PrintHitsInfo();
00814         }
00815       }
00816 
00817       // We found all possible 2D roads in one part and seg
00818       // and now it is time to construct 3D roads base on those 2D roads
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           // Create a new road object with road0 and road1
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           // Check that both 1D roads are in the same part and segment.
00834           if ( (road0->GetPart() != road1->GetPart()) 
00835               || (road0->GetSeg() != road1->GetSeg()) ) {
00836             continue;
00837           }
00838 
00839           MucRec3DRoad *road3D = new MucRec3DRoad(road0, road1);
00840 
00841           // What are our criteria for accepting a road as valid?
00842           // Check them here.
00843           // 1. difference of total number clusters in road0 and in road1
00844           // 2. difference of last gaps in road0 and road1
00845           // 3. InterSection OK or not (if X Y clusters in each gap are
00846           //                            in same panel or overlap region) 
00847           // 4. Reference position check
00848           // 5. Chisquare check
00849           // 6. Last Gap check
00850           // 7. Duplicate road check
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             //cout << "TotalHitsDelta = " << road3D->GetTotalHitsDelta() << endl;
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             //cout<<"judge shared hits:  road["<<i<<"]   max = "<<maxSharedHits<<endl;
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);//gap0
00931 
00932             //cout<<"slope1,0 = "<<slope1<<"  "<<slope0<<"  vx,y = "<<vx<<"  "<<vy<<endl;
00933             //cout<<"startxyz= "<<startx<<"  "<<starty<<"  "<<startz<<endl;
00934             //mom(vx,vy,1)
00935             float vz = 1;
00936             float sign = vy/fabs(vy);
00937             float signx = vx/fabs(vx);
00938             //cout<<"vxyz = "<<vx<<" "<<vy<<" "<<vz<<endl;
00939             if(road3D->GetPart() == 1){
00940               if(road3D->GetSeg() > 4){   //down segment
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               //fix me
00964 
00965               //cout<<"startxyz= "<<startx<<"  "<<starty<<"  "<<startz<<endl;
00966               //cout<<"vx,y,z = "<<vx<<"  "<<vy<<"  "<<vz<<endl;
00967               //cout<<"in road finder a endcap finded!!! -------------"<<endl;
00968             }
00969             else if(road3D->GetPart() == 2){
00970               //fix me
00971 
00972               vx *= -1;
00973               vy *= -1;
00974               vz *= -1;
00975               //cout<<"startxyz= "<<startx<<"  "<<starty<<"  "<<startz<<endl;
00976               //cout<<"vx,y,z = "<<vx<<"  "<<vy<<"  "<<vz<<endl;
00977               //cout<<"in road finder a endcap finded!!! ------------2-"<<endl;
00978 
00979             }
00980 
00981 
00982             Hep3Vector mom(vx, vy, vz);
00983 
00985             //MucTrack *aTrack = new MucTrack(road3D);
00986             //cout<<"startxyz = "<<startx<<"  "<<starty<<"  "<<startz<<endl;
00987             //cout<<"vxyz = "<<vx<<"  "<<vy<<"  "<<vz<<endl;
00988 
00989             startx /= 10; starty /= 10; startz /= 10;  //mm->cm
00990             startx -= vx/mom.mag(); starty -= vy/mom.mag(); startz -= vz/mom.mag();  // decrease a little  
00991 
00992             //cout<<"startxyz = "<<startx<<"  "<<starty<<"  "<<startz<<endl;
00993             RecMucTrack *aTrack = new RecMucTrack();
00994             aTrack->SetExtMucPos(startx, starty, startz);  //mm->cm
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             //aTrack->LineFit(1);
01002             //aTrack->ComputeTrackInfo(1);
01003 
01004             aTrack->setTrackId(trackId);
01005             aTrack->setId(muctrackId);
01006             trackId++;
01007             muctrackId++;
01008 
01009             //cout<<"find a track!!!"<<endl; 
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               //cout<<"expected Hits: "<<ihit->Part()<<"  "<<ihit->Seg()<<"  "<<ihit->Gap()<<"  "<<ihit->Strip()<<endl;
01024             }
01025 
01026             vector<int> multiHit;
01027             for(int i=0; i< attachedHits.size(); i++)
01028             {
01029               MucRecHit *ihit = attachedHits[i];
01030               //cout<<"attached Hits: "<<ihit->Part()<<"  "<<ihit->Seg()<<"  "<<ihit->Gap()<<"  "<<ihit->Strip()<<" hitmode: "<<ihit->GetHitMode()<<endl;
01031               //calc multiplicity hits;
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               //cout<<"multi hit: "<<hitnum<<" "<<multiHit[i]<<endl;
01040 
01041             }
01042 
01043             for(int i=0; i< expectedHits.size(); i++)
01044             { //calc distance between attached hits and expected hits
01045 
01046               MucRecHit *ihit = expectedHits[i];
01047               //cout<<"attached Hits: "<<ihit->Part()<<"  "<<ihit->Seg()<<"  "<<ihit->Gap()<<"  "<<ihit->Strip()<<endl;
01048 
01049               int diff = -999;
01050 
01051               for(int j=0; j< attachedHits.size(); j++){
01052                 MucRecHit *jhit = attachedHits[j];
01053 
01054                 //                if(attachedHits.size()!=distanceHits.size())cout<<"attached hits size no same as disthits!!!"<<endl;
01055 
01056                 if((jhit->Part()==ihit->Part())&&(jhit->Seg()==ihit->Seg())&&(jhit->Gap()==ihit->Gap())&&attachedHits.size()==distanceHits.size())
01057                 {    // same gap, cale distance
01058                   diff = ihit->Strip() - jhit->Strip();
01059                   //cout<<"diff = "<<diff<<endl;
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                     //cout<<"distance = "<<m_dist<<endl;
01068 
01069                     m_angle_cosmic = -999;
01070                     m_angle_updown = -999;
01071                     //m_px = -999; m_py = -999; m_pz = -999; m_theta = -999; m_phi = -999; m_theta_pipe = -999; m_phi_pipe = -999;
01072                     //m_px_mc = -999; m_py_mc = -999; m_pz_mc = -999; m_theta_mc = -999; m_phi_mc = -999;
01073                     //m_theta_mc_pipe = -999; m_phi_mc_pipe = -999;
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) {  // to calc effi of this strip
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                   //m_px = -999; m_py = -999; m_pz = -999; m_theta = -999; m_phi = -999; m_theta_pipe = -999; m_phi_pipe = -999;
01099                   //m_px_mc = -999; m_py_mc = -999; m_pz_mc = -999; m_theta_mc = -999; m_phi_mc = -999;
01100                   //m_theta_mc_pipe = -999; m_phi_mc_pipe = -999;
01101                   m_ngapwithhits = aTrack->numLayers();
01102                   m_run = eventHeader->runNumber();
01103                   m_event = eventHeader->eventNumber();
01104                   //m_tuple->write();
01105                 }
01106               }
01107               //if(diff != -999) cout<< "has hit in this layer"<<endl;
01108 
01109             }
01111             /*
01112                m_part = -999;
01113                m_seg  = -999;
01114                m_gap  = -999;
01115                m_strip= -999;
01116                m_diff = -999;
01117 
01118                m_angle_updown = -999;
01119                m_angle_cosmic =  cosmicMom.angle(aTrack->getMucMomentum());
01120                if(m_angle_cosmic>1.57) m_angle_cosmic = 3.14159265 - m_angle_cosmic;
01121                m_run = eventHeader->runNumber();
01122                m_event = eventHeader->eventNumber();
01123 
01124                double px,py,pz,p,theta,phi;
01125                px = (aTrack->getMucMomentum()).x(); 
01126                py = (aTrack->getMucMomentum()).y();
01127                pz = (aTrack->getMucMomentum()).z();
01128 
01129                Hep3Vector rotate(-px,pz,py);
01130                theta = rotate.theta();
01131                phi   = rotate.phi();          
01132 
01133 
01134                m_px = px; m_py = py; m_pz = pz;
01135                m_theta = theta; m_phi = phi; 
01136                m_theta_pipe = (aTrack->getMucMomentum()).theta();
01137                m_phi_pipe   = (aTrack->getMucMomentum()).phi();
01138             //m_px_mc = -999; m_py_mc = -999; m_pz_mc = -999; m_theta_mc = -999; m_phi_mc = -999;
01139             //m_theta_mc_pipe = -999; m_phi_mc_pipe = -999;
01141 
01142              */
01143 
01144           } 
01145           else {
01146             //cout << "Delete a 3D Road ... " << endl;
01147             delete road3D;
01148             // don't keep it if it's not a good condidate
01149           }
01150         } // for ( int iRoad1 ...
01151       } // for ( int iRoad0 .. 
01152 
01153 
01154       //for cosmic ray, to combine 2 track to 1
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           //if(fabs(m_phi_pipe*57.2958-90)>3)cout<<"px,y,z = "<<m_px<<"  "<<m_py<<" "<<m_pz<<endl;
01191 
01192           //calc proj point in y=0 plane
01193           Hep3Vector mucPos = aaTrack->getMucPos();
01194           double posx, posy, posz;
01195           posx = mucPos.x(); posy = mucPos.y(); posz = mucPos.z();
01196           //double projx, projz;
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           //cout<<"projection: "<<projx<<"  "<<projz<<endl;
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;   // this track has been used
01217           aaTrack->setTrackId(iTrack);
01218           //cout<<"atrack set index "<<iTrack<<endl;
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             //cout<<"angle is : "<<mom_a.angle(mom_b)<<endl;
01226             if(fabs(mom_a.angle(mom_b) - 3.1415926)<0.2)
01227             {
01228               bbTrack->setTrackId(iTrack);
01229               //cout<<"btrack set index "<<iTrack<<endl;
01230               //cout<<"find one cosmic ray"<<endl;
01231 
01232               //cout<<"angle = "<<cosmicMom.angle(mom_a)<<"  "<<cosmicMom.angle(mom_b)<<endl;
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               //m_tuple->write();
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             //road->PrintHitsInfo();
01278           }
01279         }
01280 
01281         m_NEventReco ++;
01282       }
01283 
01284       //delete p3DRoadC
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 }

StatusCode MucRecRoadFinder::finalize (  ) 

Definition at line 1302 of file MucRecRoadFinder.cxx.

References Bes_Common::INFO, and msgSvc().

01303 {
01304   MsgStream log(msgSvc(), name());
01305   log << MSG::INFO << "in finalize()" << endreq;
01306 
01307   //cout << m_NHitsLostTotal << " of " << m_NHitsTotal << " total hits" << endl; 
01308   //for(int i = 0; i < 20; i++) cout << "lost " << i << " hits event " << m_NHitsLost[i] << endl;
01309   //for(int i = 0; i < 9; i++) cout << "lost on gap " << i << " is " << m_NHitsLostInGap[i] << endl;
01310 
01311   return StatusCode::SUCCESS;
01312 }

StatusCode MucRecRoadFinder::initialize (  ) 

Definition at line 66 of file MucRecRoadFinder.cxx.

References aMucRecHitContainer, IMucGeomSvc::Dump(), calibUtil::ERROR, genRecEmupikp::i, Bes_Common::INFO, m_angle_cosmic, m_angle_updown, m_diff, m_dist, m_emcDown, m_emcUp, m_event, m_filter_event, m_filter_filename, m_gap, m_maxhit, m_mucDown, m_mucUp, m_multihit, m_NEvent, m_NEventReco, m_NEventWithHit, m_ngapwithhits, m_nhit, m_NHitsLost, m_NHitsLostInGap, m_NHitsLostTotal, m_NHitsTotal, m_NtOutput, 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_seg, m_strip, m_theta, m_theta_mc, m_theta_mc_pipe, m_theta_pipe, m_tuple, msgSvc(), and ntupleSvc().

00067 {
00068   MsgStream log(msgSvc(), name());
00069   log << MSG::INFO << "in initialize()" << endreq;
00070   
00071   m_NHitsTotal = 0;
00072   m_NHitsLostTotal = 0;
00073   for(int i = 0; i < 20; i++) m_NHitsLost.push_back(0);
00074   for(int i = 0; i < 10; i++) m_NHitsLostInGap.push_back(0);
00075 
00076   m_NEvent        = 0;
00077   m_NEventWithHit = 0;
00078   m_NEventReco    = 0;
00079 
00080   IMucGeomSvc* mucGeomSvc;
00081   StatusCode sc = service("MucGeomSvc", mucGeomSvc);
00082   if (sc ==  StatusCode::SUCCESS) {              
00083     mucGeomSvc->Dump();
00084     //cout<<"1st wire id:"<<mucGeomSvc->Wire(0)->Id()<<endl;
00085     //cout<<"2nd wire lyr id:"<<mucGeomSvc->Wire(1)->Lyr()->Id()<<endl;
00086     //cout<<"6860th wire lyr id:"<<mucGeomSvc->Wire(6859)->Lyr()->Id()<<endl; 
00087   } else {
00088     return StatusCode::FAILURE;
00089   }
00090 
00091   aMucRecHitContainer = 0;
00092 
00093   if(m_NtOutput>=1){ 
00094   NTuplePtr nt1(ntupleSvc(), "FILE401/T");
00095   
00096   if ( nt1 ) { m_tuple = nt1;}
00097   else { 
00098     //    m_tuple = ntupleSvc()->book ("MyTuples/1", CLID_RowWiseTuple, "MdcTrkRecon N-Tuple");
00099     m_tuple = ntupleSvc()->book ("FILE401/T", CLID_RowWiseTuple, "MucTrkRecon N-Tuple");
00100     if ( m_tuple )    {
00101       sc = m_tuple->addItem ("part",  m_part);
00102       sc = m_tuple->addItem ("seg",   m_seg);
00103       sc = m_tuple->addItem ("gap",   m_gap);
00104       sc = m_tuple->addItem ("strip", m_strip);
00105       sc = m_tuple->addItem ("diff",  m_diff);
00106       sc = m_tuple->addItem ("dist",  m_dist);
00107       sc = m_tuple->addItem ("run",   m_run);
00108       sc = m_tuple->addItem ("event", m_event);
00109       sc = m_tuple->addItem ("ngap",  m_ngapwithhits);
00110       sc = m_tuple->addItem ("nhit",  m_nhit);
00111       sc = m_tuple->addItem ("maxhit",  m_maxhit);
00112       sc = m_tuple->addItem ("multihit",m_multihit);
00113       sc = m_tuple->addItem ("angleCosmic",m_angle_cosmic);
00114       sc = m_tuple->addItem ("angleUpdown",m_angle_updown);  
00115       sc = m_tuple->addItem ("px",m_px);
00116       sc = m_tuple->addItem ("py",m_py);
00117       sc = m_tuple->addItem ("pz",m_pz);
00118       sc = m_tuple->addItem ("theta",m_theta);
00119       sc = m_tuple->addItem ("phi",m_phi);
00120       sc = m_tuple->addItem ("theta_pipe",m_theta_pipe);
00121       sc = m_tuple->addItem ("phi_pipe",m_phi_pipe);
00122       sc = m_tuple->addItem ("pxmc",m_px_mc);
00123       sc = m_tuple->addItem ("pymc",m_py_mc);
00124       sc = m_tuple->addItem ("pzmc",m_pz_mc);
00125       sc = m_tuple->addItem ("thetamc",m_theta_mc);
00126       sc = m_tuple->addItem ("phimc",m_phi_mc);
00127       sc = m_tuple->addItem ("thetamc_pipe",m_theta_mc_pipe);
00128       sc = m_tuple->addItem ("phimc_pipe",m_phi_mc_pipe);
00129       sc = m_tuple->addItem ("emcUp",m_emcUp);
00130       sc = m_tuple->addItem ("emcDown",m_emcDown);
00131       sc = m_tuple->addItem ("mucUp",m_mucUp);
00132       sc = m_tuple->addItem ("mucDown",m_mucDown);
00133       sc = m_tuple->addItem ("projx",m_projx);
00134       sc = m_tuple->addItem ("projz",m_projz);
00135    
00136     }
00137     else    {   // did not manage to book the N tuple....
00138       log << MSG::ERROR << "    Cannot book N-tuple:" << long(m_tuple) << endmsg;
00139       //return StatusCode::FAILURE;
00140     }
00141   }
00142   }
00143 
00144   // load the filter event
00145   if ( m_filter_filename == "") {
00146    m_filter_filename =getenv("MUCRECALGROOT");
00147    m_filter_filename +="/share/filter.txt";
00148    }
00149 
00150   if (m_filter_filename.size()) {
00151     std::ifstream infile(m_filter_filename.c_str());
00152 
00153     while (!infile.eof()) {
00154       FilterEvent filterevt;
00155       infile >> filterevt.bossver >> filterevt.runid >> filterevt.eventid;
00156       if ((!infile.good()) || (infile.eof())) {
00157         break;
00158       }
00159 
00160       m_filter_event.push_back(filterevt);
00161 //      cout << filterevt.bossver << " " 
00162 //           << filterevt.runid << " "
00163 //           << filterevt.eventid << std::endl;
00164     }
00165   }
00166   
00167   return StatusCode::SUCCESS;
00168 }

void MucRecRoadFinder::TrackFinding ( RecMucTrack aTrack  ) 

Definition at line 1315 of file MucRecRoadFinder.cxx.

References aMucRecHitContainer, RecMucTrack::AttachHit(), RecMucTrack::AttachInsct(), RecMucTrack::CalculateInsct(), RecMucTrack::ComputeTrackInfo(), RecMucTrack::CorrectDir(), RecMucTrack::CorrectPos(), Bes_Common::DEBUG, RecMucTrack::GetCurrentDir(), RecMucTrack::GetCurrentInsct(), RecMucTrack::GetCurrentPos(), MucGeoGeneral::GetGap(), MucRecHitContainer::GetGapHitCount(), MucID::getGapNum(), MucRecHitContainer::GetHit(), RecMucTrack::GetHitDistance(), MucRecHit::GetHitMode(), MucID::getPartNum(), RecMucTrack::GetRecMode(), MucID::getSegNum(), RecMucTrack::GetTotalHits(), DstMucTrack::id(), Bes_Common::INFO, MucGeoGeneral::Instance(), kDeltaSeg, kNSegSearch, kPartSeq, kWindowSize, RecMucTrack::LineFit(), m_MsOutput, m_NHitsLostInGap, m_onlyseedfit, msgSvc(), MucGeoGap::Orient(), RecMucTrack::Project(), RecMucTrack::SetCurrentInsct(), MucRecHit::SetHitMode(), MucRecHit::Strip(), DstMucTrack::trackId(), and Bes_Common::WARNING.

Referenced by execute().

01316 {
01317   MsgStream log(msgSvc(), name());
01318 
01319   Hep3Vector currentPos = aTrack->GetCurrentPos();
01320   Hep3Vector currentDir = aTrack->GetCurrentDir();
01321   //   if(currentPos.mag() < kMinor) {
01322   //     log << MSG::WARNING << "No MUC intersection in track " << endreq;
01323   //     continue;
01324   //   }
01325 
01326   int firstHitFound[2] = { 0,  0}; // Has the fist position in this orient determined? if so, could CorrectDirection()
01327   int firstHitGap[2]   = {-1, -1}; // When the fist position in this orient determined, the gap it is on
01328   for(int partSeq = 0; partSeq < (int)MucID::getPartNum(); partSeq++) {
01329     int iPart = kPartSeq[partSeq];
01330     for(int iGap = 0; iGap < (int)MucID::getGapNum(iPart); iGap++) {
01331       int seg = -1;
01332       int orient = MucGeoGeneral::Instance()->GetGap(iPart, 0, iGap)->Orient();;
01333 
01334       float xInsct, yInsct, zInsct;
01335       aTrack->Project(iPart, iGap, xInsct, yInsct, zInsct, seg);
01336       if(m_MsOutput) cout<<"part "<<iPart<<" gap "<<iGap<<" x "<<xInsct<<" y "<<yInsct<<" z "<<zInsct<<" seg "<<seg<<endl; 
01337 
01338       if(seg == -1) {
01339         //log << MSG::DEBUG << "no intersection with part " << iPart
01340         //  << " gap " << iGap << " in this track !" << endl;
01341         continue;
01342       }
01343 
01344       aTrack->SetCurrentInsct(xInsct, yInsct, zInsct);
01345 
01346       for(int iDeltaSeg = 0; iDeltaSeg < kNSegSearch; iDeltaSeg++) {
01347         int iSeg = seg + kDeltaSeg[iDeltaSeg];  // also find on neighbor seg;
01348         if(iSeg < 0) iSeg += MucID::getSegNum(iPart);
01349         if(iSeg > (int)MucID::getSegNum(iPart) - 1) iSeg -= MucID::getSegNum(iPart);
01350 
01351         float window = kWindowSize[iPart][iGap];
01352 
01353         for (int iHit = 0; iHit < aMucRecHitContainer->GetGapHitCount(iPart, iSeg, iGap); iHit++) {
01354           log << MSG::DEBUG << "iSeg " << iSeg << " iHit " << iHit << endreq;
01355           MucRecHit* pHit = aMucRecHitContainer->GetHit(iPart, iSeg, iGap, iHit);
01356           //cout<< "strip " << pHit->Strip() << " center " << pHit->GetCenterPos() << endl;
01357 
01358           if (!pHit) {
01359             log << MSG::WARNING << "MucRecTrkExt: null pointer to pHit" << endreq;
01360             continue;
01361           }
01362           else {
01363 
01364             //this algo use with ext and this hit has been used before;                   
01365             if(pHit->GetHitMode() != -1 && pHit->GetHitMode() != 3) continue;
01366 
01367             // Get the displacement of hit pHit to intersection
01368             float dX  = aTrack->GetHitDistance(pHit);
01369             log << MSG::DEBUG << "distance = " << setw(8) << dX << " size " << setw(4) << window << endreq;
01370 
01371             //cout<<"dX= "<<dX<<"  window="<<window<<endl;
01372 
01373             if (dX < window) {
01374               // Attach the hit if it exists
01375               //cout << "meet window " << pHit->GetSoftID() << endl;
01376               //****************if this if emc track, we abondon used hit in mdc*******************
01377               //if(m_emcrec == 1 )
01378               if(aTrack->GetRecMode() == 0) {
01379                 pHit->SetHitMode(1); //mdc ext
01380                 aTrack->AttachHit(pHit);
01381                 //cout<<"in MucRecTrkExt: trackmode==0 so mdc ext "<<iPart<<" "<<iSeg<<" "<<iGap<<" "<<iHit<<endl;
01382               }
01383               if(aTrack->GetRecMode() == 1) {
01384                 //cout<<"in MucRecTrkExt: HitMode = "<<pHit->GetHitMode()<<" part:  "<<iPart<<" "<<iSeg<<" "<<iGap<<" "<<iHit<<endl;
01385                 if(pHit->GetHitMode()!=1) {
01386                   aTrack->AttachHit(pHit);  //this hit has not been used by mdc ext
01387                   pHit->SetHitMode(2); //emc ext
01388                 }
01389               }
01390               if(aTrack->GetRecMode() == 2) {
01391                 aTrack->AttachHit(pHit);  //this hit has not been used by mdc ext
01392                 //pHit->SetHitMode(2); //emc ext
01393               }
01394               if(aTrack->GetRecMode() == 3) { 
01395                 if(pHit->GetHitMode()!=1) {     
01396                   aTrack->AttachHit(pHit);  //this hit has not been used by mdc ext
01397                   pHit->SetHitMode(3); //emc ext
01398                 }
01399               }
01400 
01401               if (firstHitGap[orient] == -1) firstHitGap[orient] = iGap;
01402               firstHitFound[orient] = 1;
01403               //cout << "You could correct directon in orient " << orient << endl;
01404 
01405               //cout<< " part " << iPart << " seg " << iSeg << " gap " << iGap 
01406               //          << " strip " << setw(2) << pHit->Strip() << " attatched" << endl;
01407 
01408               log << MSG::DEBUG << " part " << iPart << " seg " << iSeg << " gap " << iGap 
01409                 << " strip " << setw(2) << pHit->Strip() << " attatched" << endreq;
01410               log << MSG::DEBUG << "current total hits " << aTrack->GetTotalHits() << endreq;
01411             }
01412             else {
01413               m_NHitsLostInGap[iGap]++;
01414               //log << MSG::DEBUG cout << " part " << iPart << " seg " << iSeg << " gap " << iGap 
01415               //       << " strip " << setw(2) << pHit->GetSoftID().GetStrip() 
01416               //     << " not attached !" << endreq;
01417             }
01418           }
01419         }
01420         aTrack->CalculateInsct(iPart, iSeg, iGap);
01421       }
01422 
01423       if(m_onlyseedfit == 0) {//
01424         // When correct dir in the orient is permitted and this gap is not the gap first hit locates. 
01425         if (firstHitFound[orient] && firstHitGap[orient] != iGap) aTrack->CorrectDir();
01426         aTrack->CorrectPos();
01427         //cout << "Current Intersection " << aTrack->GetCurrentInsct() << endl;
01428         //cout << "Current Direction    " << aTrack->GetCurrentDir() << endl;
01429         aTrack->AttachInsct(aTrack->GetCurrentInsct());
01430       }// else; not correct pos & dir.
01431 
01432     }
01433   }
01434   aTrack->LineFit(1);
01435   aTrack->ComputeTrackInfo(1);
01436   log << MSG::INFO << "Fit track done! trackIndex: " << aTrack->trackId() << ", mucId: " << aTrack->id() << ", RecMode: " << aTrack->GetRecMode() << endreq; 
01437   //cout<<" in TrackFinding: depth= "<<aTrack->GetDepth()<<endl;
01438 }


Member Data Documentation

MucRecHitContainer* MucRecRoadFinder::aMucRecHitContainer [private]

Definition at line 66 of file MucRecRoadFinder.h.

Referenced by execute(), initialize(), and TrackFinding().

NTuple::Item<double> MucRecRoadFinder::m_angle_cosmic [private]

Definition at line 82 of file MucRecRoadFinder.h.

Referenced by execute(), and initialize().

NTuple::Item<double> MucRecRoadFinder::m_angle_updown [private]

Definition at line 83 of file MucRecRoadFinder.h.

Referenced by execute(), and initialize().

std::string MucRecRoadFinder::m_configFile [private]

Definition at line 50 of file MucRecRoadFinder.h.

Referenced by MucRecRoadFinder().

NTuple::Item<double> MucRecRoadFinder::m_diff [private]

Definition at line 73 of file MucRecRoadFinder.h.

Referenced by execute(), and initialize().

NTuple::Item<double> MucRecRoadFinder::m_dist [private]

Definition at line 74 of file MucRecRoadFinder.h.

Referenced by execute(), and initialize().

NTuple::Item<double> MucRecRoadFinder::m_emcDown [private]

Definition at line 99 of file MucRecRoadFinder.h.

Referenced by execute(), and initialize().

NTuple::Item<double> MucRecRoadFinder::m_emcUp [private]

Definition at line 98 of file MucRecRoadFinder.h.

Referenced by execute(), and initialize().

NTuple::Item<double> MucRecRoadFinder::m_event [private]

Definition at line 77 of file MucRecRoadFinder.h.

Referenced by execute(), and initialize().

std::vector<FilterEvent> MucRecRoadFinder::m_filter_event [private]

Definition at line 65 of file MucRecRoadFinder.h.

Referenced by execute(), and initialize().

std::string MucRecRoadFinder::m_filter_filename [private]

Definition at line 59 of file MucRecRoadFinder.h.

Referenced by initialize(), and MucRecRoadFinder().

int MucRecRoadFinder::m_fittingMethod [private]

Definition at line 49 of file MucRecRoadFinder.h.

Referenced by MucRecRoadFinder().

NTuple::Item<double> MucRecRoadFinder::m_gap [private]

Definition at line 71 of file MucRecRoadFinder.h.

Referenced by execute(), and initialize().

NTuple::Item<double> MucRecRoadFinder::m_maxhit [private]

Definition at line 80 of file MucRecRoadFinder.h.

Referenced by execute(), and initialize().

int MucRecRoadFinder::m_maxHitsRec [private]

Definition at line 54 of file MucRecRoadFinder.h.

Referenced by execute(), and MucRecRoadFinder().

int MucRecRoadFinder::m_mccosmic [private]

Definition at line 51 of file MucRecRoadFinder.h.

Referenced by execute(), and MucRecRoadFinder().

bool MucRecRoadFinder::m_MsOutput [private]

Definition at line 57 of file MucRecRoadFinder.h.

Referenced by execute(), MucRecRoadFinder(), and TrackFinding().

NTuple::Item<double> MucRecRoadFinder::m_mucDown [private]

Definition at line 101 of file MucRecRoadFinder.h.

Referenced by execute(), and initialize().

NTuple::Item<double> MucRecRoadFinder::m_mucUp [private]

Definition at line 100 of file MucRecRoadFinder.h.

Referenced by execute(), and initialize().

NTuple::Item<double> MucRecRoadFinder::m_multihit [private]

Definition at line 81 of file MucRecRoadFinder.h.

Referenced by execute(), and initialize().

int MucRecRoadFinder::m_NEvent [private]

Definition at line 41 of file MucRecRoadFinder.h.

Referenced by execute(), and initialize().

int MucRecRoadFinder::m_NEventReco [private]

Definition at line 43 of file MucRecRoadFinder.h.

Referenced by execute(), and initialize().

int MucRecRoadFinder::m_NEventWithHit [private]

Definition at line 42 of file MucRecRoadFinder.h.

Referenced by initialize().

NTuple::Item<double> MucRecRoadFinder::m_ngapwithhits [private]

Definition at line 78 of file MucRecRoadFinder.h.

Referenced by execute(), and initialize().

NTuple::Item<double> MucRecRoadFinder::m_nhit [private]

Definition at line 79 of file MucRecRoadFinder.h.

Referenced by execute(), and initialize().

std::vector<int> MucRecRoadFinder::m_NHitsLost [private]

Definition at line 46 of file MucRecRoadFinder.h.

Referenced by initialize().

std::vector<int> MucRecRoadFinder::m_NHitsLostInGap [private]

Definition at line 47 of file MucRecRoadFinder.h.

Referenced by initialize(), and TrackFinding().

int MucRecRoadFinder::m_NHitsLostTotal [private]

Definition at line 44 of file MucRecRoadFinder.h.

Referenced by initialize().

int MucRecRoadFinder::m_NHitsTotal [private]

Definition at line 45 of file MucRecRoadFinder.h.

Referenced by initialize().

int MucRecRoadFinder::m_NtOutput [private]

Definition at line 52 of file MucRecRoadFinder.h.

Referenced by execute(), initialize(), and MucRecRoadFinder().

int MucRecRoadFinder::m_onlyseedfit [private]

Definition at line 53 of file MucRecRoadFinder.h.

Referenced by execute(), MucRecRoadFinder(), and TrackFinding().

NTuple::Item<double> MucRecRoadFinder::m_part [private]

Definition at line 69 of file MucRecRoadFinder.h.

Referenced by execute(), and initialize().

NTuple::Item<double> MucRecRoadFinder::m_phi [private]

Definition at line 88 of file MucRecRoadFinder.h.

Referenced by execute(), and initialize().

NTuple::Item<double> MucRecRoadFinder::m_phi_mc [private]

Definition at line 95 of file MucRecRoadFinder.h.

Referenced by execute(), and initialize().

NTuple::Item<double> MucRecRoadFinder::m_phi_mc_pipe [private]

Definition at line 97 of file MucRecRoadFinder.h.

Referenced by execute(), and initialize().

NTuple::Item<double> MucRecRoadFinder::m_phi_pipe [private]

Definition at line 90 of file MucRecRoadFinder.h.

Referenced by execute(), and initialize().

NTuple::Item<double> MucRecRoadFinder::m_projx [private]

Definition at line 102 of file MucRecRoadFinder.h.

Referenced by execute(), and initialize().

NTuple::Item<double> MucRecRoadFinder::m_projz [private]

Definition at line 103 of file MucRecRoadFinder.h.

Referenced by execute(), and initialize().

NTuple::Item<double> MucRecRoadFinder::m_px [private]

Definition at line 84 of file MucRecRoadFinder.h.

Referenced by execute(), and initialize().

NTuple::Item<double> MucRecRoadFinder::m_px_mc [private]

Definition at line 91 of file MucRecRoadFinder.h.

Referenced by execute(), and initialize().

NTuple::Item<double> MucRecRoadFinder::m_py [private]

Definition at line 85 of file MucRecRoadFinder.h.

Referenced by execute(), and initialize().

NTuple::Item<double> MucRecRoadFinder::m_py_mc [private]

Definition at line 92 of file MucRecRoadFinder.h.

Referenced by execute(), and initialize().

NTuple::Item<double> MucRecRoadFinder::m_pz [private]

Definition at line 86 of file MucRecRoadFinder.h.

Referenced by execute(), and initialize().

NTuple::Item<double> MucRecRoadFinder::m_pz_mc [private]

Definition at line 93 of file MucRecRoadFinder.h.

Referenced by execute(), and initialize().

NTuple::Item<double> MucRecRoadFinder::m_run [private]

Definition at line 76 of file MucRecRoadFinder.h.

Referenced by execute(), and initialize().

int MucRecRoadFinder::m_seedtype [private]

Definition at line 56 of file MucRecRoadFinder.h.

Referenced by execute(), and MucRecRoadFinder().

NTuple::Item<double> MucRecRoadFinder::m_seg [private]

Definition at line 70 of file MucRecRoadFinder.h.

Referenced by execute(), and initialize().

NTuple::Item<double> MucRecRoadFinder::m_strip [private]

Definition at line 72 of file MucRecRoadFinder.h.

Referenced by execute(), and initialize().

NTuple::Item<double> MucRecRoadFinder::m_theta [private]

Definition at line 87 of file MucRecRoadFinder.h.

Referenced by execute(), and initialize().

NTuple::Item<double> MucRecRoadFinder::m_theta_mc [private]

Definition at line 94 of file MucRecRoadFinder.h.

Referenced by execute(), and initialize().

NTuple::Item<double> MucRecRoadFinder::m_theta_mc_pipe [private]

Definition at line 96 of file MucRecRoadFinder.h.

Referenced by execute(), and initialize().

NTuple::Item<double> MucRecRoadFinder::m_theta_pipe [private]

Definition at line 89 of file MucRecRoadFinder.h.

Referenced by execute(), and initialize().

NTuple::Tuple* MucRecRoadFinder::m_tuple [private]

Definition at line 68 of file MucRecRoadFinder.h.

Referenced by execute(), and initialize().

int MucRecRoadFinder::m_united [private]

Definition at line 55 of file MucRecRoadFinder.h.

Referenced by execute(), and MucRecRoadFinder().


Generated on Tue Nov 29 23:20:30 2016 for BOSS_7.0.2 by  doxygen 1.4.7