Main Page | Namespace List | Class Hierarchy | Alphabetical List | Class List | Directories | File List | Namespace Members | Class Members | File Members | Related Pages

MucRecTrkExt Class Reference

#include <MucRecTrkExt.h>

List of all members.

Public Member Functions

StatusCode execute ()
StatusCode execute ()
StatusCode finalize ()
StatusCode finalize ()
float getWindowSize (Hep3Vector mom, int part, int seg, int gap)
float getWindowSize (Hep3Vector mom, int part, int seg, int gap)
StatusCode initialize ()
StatusCode initialize ()
 MucRecTrkExt (const std::string &name, ISvcLocator *pSvcLocator)
 MucRecTrkExt (const std::string &name, ISvcLocator *pSvcLocator)
void TrackFinding (RecMucTrack *aTrack)
void TrackFinding (RecMucTrack *aTrack)

Private Attributes

bool m_Blind
NTuple::Item< double > m_Chi2
NTuple::Item< double > m_Chi2
int m_CompareWithMcHit
std::string m_configFile
NTuple::Item< double > m_depth
NTuple::Item< double > m_depth
NTuple::Item< double > m_diff
NTuple::Item< double > m_diff
NTuple::Item< double > m_Dist_x
NTuple::Item< double > m_Dist_x
NTuple::Item< double > m_Dist_y
NTuple::Item< double > m_Dist_y
NTuple::Item< double > m_Dist_z
NTuple::Item< double > m_Dist_z
NTuple::Item< double > m_distance
NTuple::Item< double > m_Distance
NTuple::Item< double > m_distance
NTuple::Item< double > m_Distance
int m_emcrec
int m_EmcShowerSeedMode
NTuple::Item< double > m_emctrack
NTuple::Item< double > m_emctrack
NTuple::Item< double > m_event
NTuple::Item< double > m_event
int m_ExtTrackSeedMode
int m_fittingMethod
NTuple::Item< double > m_gap
NTuple::Item< double > m_gap
NTuple::Item< double > m_MaxHits
NTuple::Item< double > m_MaxHits
bool m_MsOutput
NTuple::Item< double > m_muc_digi
NTuple::Item< double > m_muc_digi
int m_MucHitSeedMode
MucRecHitContainer * m_MucRecHitContainer
MucRecHitContainer * m_MucRecHitContainer
int m_NDigisTotal
int m_NHitsFoundTotal
std::vector< int > m_NHitsLost
std::vector< int > m_NHitsLost
int m_NHitsLostByExtTotal
int m_NHitsLostByMdcTotal
std::vector< int > m_NHitsLostInGap
std::vector< int > m_NHitsLostInGap
int m_NHitsLostTotal
int m_NHitsMisFoundTotal
int m_NHitsTotal
int m_NtOutput
int m_NTracksLostByExtTotal
int m_NTracksLostByMdcTotal
int m_NTracksMdcGhostTotal
int m_NTracksRecoTotal
int m_NTracksTotal
NTuple::Item< double > m_part
NTuple::Item< double > m_part
HepPDT::ParticleDataTable * m_particleTable
HepPDT::ParticleDataTable * m_particleTable
NTuple::Item< double > m_posx
NTuple::Item< double > m_posx
NTuple::Item< double > m_posx_ext
NTuple::Item< double > m_posx_ext
NTuple::Item< double > m_posx_sigma
NTuple::Item< double > m_posx_sigma
NTuple::Item< double > m_posy
NTuple::Item< double > m_posy
NTuple::Item< double > m_posy_ext
NTuple::Item< double > m_posy_ext
NTuple::Item< double > m_posy_sigma
NTuple::Item< double > m_posy_sigma
NTuple::Item< double > m_posz
NTuple::Item< double > m_posz
NTuple::Item< double > m_posz_ext
NTuple::Item< double > m_posz_ext
NTuple::Item< double > m_posz_sigma
NTuple::Item< double > m_posz_sigma
NTuple::Item< double > m_px_mc
NTuple::Item< double > m_px_mc
NTuple::Item< double > m_py_mc
NTuple::Item< double > m_py_mc
NTuple::Item< double > m_pz_mc
NTuple::Item< double > m_pz_mc
NTuple::Item< double > m_run
NTuple::Item< double > m_run
NTuple::Item< double > m_seg
NTuple::Item< double > m_seg
NTuple::Item< double > m_strip
NTuple::Item< double > m_strip
long m_totalEvent
std::vector< int > m_TrackLostByExt
std::vector< int > m_TrackLostByExt
std::vector< int > m_TrackLostByMdc
std::vector< int > m_TrackLostByMdc
std::vector< int > m_TrackLostHit
std::vector< int > m_TrackLostHit
std::vector< int > m_TrackLostHitNumber
std::vector< int > m_TrackLostHitNumber
NTuple::Tuple * m_tuple
NTuple::Tuple * m_tuple

Constructor & Destructor Documentation

MucRecTrkExt::MucRecTrkExt (  const std::string &  name,
ISvcLocator *  pSvcLocator
) 
  

00048                                                                           :
00049   Algorithm(name, pSvcLocator)
00050 {
00051   // Declare the properties  
00052   declareProperty("ExtTrackSeedMode", m_ExtTrackSeedMode = 1); // 1: My ext from Mc track, 2: from Ext track, 3: from MUC for cosmic ray
00053   declareProperty("CompareWithMcHit", m_CompareWithMcHit = 0); 
00054   declareProperty("FittingMethod",    m_fittingMethod = 1);
00055   declareProperty("EmcShowerSeedMode",m_EmcShowerSeedMode = 0); // 0: Not use EmcShower as seed, 1: use...
00056   declareProperty("MucHitSeedMode",   m_MucHitSeedMode = 0); // 0: Not use MucHit as seed,1 use... 
00057   declareProperty("ConfigFile",       m_configFile = "MucConfig.xml");
00058   declareProperty("Blind",            m_Blind = false);
00059   declareProperty("NtOutput",         m_NtOutput = 0); // NTuple save to root or not
00060   declareProperty("MsOutput",         m_MsOutput = false); // Debug message cout or not
00061   
00062 }

MucRecTrkExt::MucRecTrkExt (  const std::string &  name,
ISvcLocator *  pSvcLocator
) 
 

Member Function Documentation

StatusCode MucRecTrkExt::execute (   ) 
 

StatusCode MucRecTrkExt::execute (   ) 
  

00164                                  {
00165 
00166   MsgStream log(msgSvc(), name());
00167   log << MSG::INFO << "in execute()" << endreq;
00168   
00169   
00170   StatusCode sc = StatusCode::SUCCESS;
00171 
00172   // Part 1: Get the event header, print out event and run number
00173   int event, run;
00174 
00175   SmartDataPtr<Event::EventHeader> eventHeader(eventSvc(),"/Event/EventHeader");
00176   if (!eventHeader) {
00177     log << MSG::FATAL << "Could not find Event Header" << endreq;
00178     //return( StatusCode::FAILURE);
00179   }
00180   m_totalEvent ++;
00181   log << MSG::INFO << "Event: "<<m_totalEvent<<"\tcurrent run: "<<eventHeader->runNumber()<<"\tcurrent event: "<< eventHeader->eventNumber()<<endreq;
00182 
00183   //Part 2: Retrieve MC truth 
00184   if(m_CompareWithMcHit==1){
00185     SmartDataPtr<McParticleCol> mcParticleCol(eventSvc(),"/Event/MC/McParticleCol");
00186     if (!mcParticleCol) {
00187       log << MSG::FATAL << "Could not find McParticle" << endreq;
00188       //return( StatusCode::FAILURE);
00189     }
00190 
00191     McParticleCol::iterator iter_mc = mcParticleCol->begin();
00192     //do not loop; just get first particle info
00193 
00194     //if(!(*iter_mc)->primaryParticle()) continue;
00195     int pid = (*iter_mc)->particleProperty();
00196     int charge = 0;
00197     double mass = -1.0; 
00198 
00199     if( pid >0 ) 
00200     { 
00201       if(m_particleTable->particle( pid )){
00202         charge = (int)m_particleTable->particle( pid )->charge();
00203         mass = m_particleTable->particle( pid )->mass();
00204       }
00205       else log << MSG::ERROR << "no this particle id in particle table, please check data" <<endreq;
00206     } else if ( pid <0 ) {
00207       if(m_particleTable->particle( -pid )){
00208         charge = (int)m_particleTable->particle( -pid )->charge();
00209         charge *= -1;
00210         mass = m_particleTable->particle( -pid )->mass();
00211       }
00212       else log << MSG::ERROR << "no this particle id in particle table, please check data" <<endreq;
00213     } else {
00214       log << MSG::WARNING << "wrong particle id, please check data" <<endreq;
00215     }
00216 
00217     //  if(!charge) {
00218     //    log << MSG::WARNING
00219     //        << "neutral particle charge = 0!!! ...just skip it ! " << endreq;
00220     //    continue;
00221     //  }
00222 
00223     HepLorentzVector initialMomentum = (*iter_mc)->initialFourMomentum();
00224     HepLorentzVector initialPos = (*iter_mc)->initialPosition();
00225     if(m_NtOutput>=1){
00226       m_px_mc = initialMomentum.px();
00227       m_py_mc = initialMomentum.py();
00228       m_pz_mc = initialMomentum.pz();
00229     }
00230     //cout<<"mc mom: "<<m_px_mc<<endl;
00231 
00232     log << MSG::INFO << " particleId = " << (*iter_mc)->particleProperty() << endreq;   
00233   }
00234 
00235   //Part 6: Retrieve MUC digi 
00236   int digiId = 0;
00237   SmartDataPtr<MucDigiCol> mucDigiCol(eventSvc(),"/Event/Digi/MucDigiCol");
00238   if (!mucDigiCol) {
00239     log << MSG::FATAL << "Could not find MUC digi" << endreq;
00240     return( StatusCode::FAILURE);
00241   }
00242   MucDigiCol::iterator iter4 = mucDigiCol->begin();
00243   for (;iter4 != mucDigiCol->end(); iter4++, digiId++) {
00244   }
00245   log << MSG::INFO << "Total MUC digis:\t" << digiId << endreq;
00246   if( digiId == 0 ) return ( StatusCode::SUCCESS);
00247 
00248   // Retrieve MdcMcHit
00249   /*
00250      SmartDataPtr<Event::MdcMcHitCol> aMdcMcHitCol(eventSvc(),"/Event/MC/MdcMcHitCol");
00251      if (!aMdcMcHitCol) {
00252      log << MSG::WARNING << "Could not find MdcMcHitCol" << endreq;
00253   //return( StatusCode::FAILURE);
00254   }
00255   //log << MSG::DEBUG << "MdcMcHitCol contains " << aMdcMcHitCol->size() << " Hits " << endreq;
00256 
00257   // Retrieve TofMcHit 
00258   SmartDataPtr<Event::TofMcHitCol> aTofMcHitCol(eventSvc(),"/Event/MC/TofMcHitCol");
00259   if (!aTofMcHitCol) {
00260   log << MSG::WARNING << "Could not find TofMcHitCol" << endreq;
00261   //return( StatusCode::FAILURE);
00262   }
00263   //log << MSG::DEBUG << "TofMcHitCol contains " << aTofMcHitCol->size() << " Hits " << endreq;
00264 
00265   // Retrieve EmcMcHit 
00266   SmartDataPtr<Event::EmcMcHitCol> aEmcMcHitCol(eventSvc(),"/Event/MC/EmcMcHitCol");
00267   if (!aEmcMcHitCol) {
00268   log << MSG::WARNING << "Could not find EmcMcHitCol" << endreq;
00269   //return( StatusCode::FAILURE);
00270   }
00271   //log << MSG::DEBUG << "EmcMcHitCol contains " << aEmcMcHitCol->size() << " Hits " << endreq;
00272    */
00273 
00274   Identifier mucID;
00275 
00276   //McPartToMucHitTab assocMcMuc;
00277   //assocMcMuc.init();
00278 
00279   if (m_CompareWithMcHit) {
00280     McPartToMucHitTab assocMcMuc;
00281     assocMcMuc.init();
00282 
00283     SmartDataPtr<McParticleCol> mcParticleCol(eventSvc(),"/Event/MC/McParticleCol");
00284     if (!mcParticleCol) {
00285       log << MSG::FATAL << "Could not find McParticle" << endreq;
00286       //return( StatusCode::FAILURE);
00287     }
00288     McParticleCol::iterator iter_mc = mcParticleCol->begin();
00289 
00290     // Retrieve MucMcHit 
00291     SmartDataPtr<Event::MucMcHitCol> aMucMcHitCol(eventSvc(),"/Event/MC/MucMcHitCol");
00292     if (!aMucMcHitCol) {
00293       log << MSG::WARNING << "Could not find MucMcHitCol" << endreq;
00294       //return( StatusCode::FAILURE);
00295     }
00296 
00297     log << MSG::DEBUG << "MucMcHitCol contains " << aMucMcHitCol->size() << " Hits " << endreq;
00298     MucMcHitCol::iterator iter_MucMcHit = aMucMcHitCol->begin();
00299     for (; iter_MucMcHit != aMucMcHitCol->end(); iter_MucMcHit++) {
00300       mucID = (*iter_MucMcHit)->identify();
00301 
00302       log << MSG::DEBUG << " MucMcHit " << " : " 
00303         << " part " << MucID::barrel_ec(mucID)
00304         << " seg "  << MucID::segment(mucID)
00305         << " gap "  << MucID::layer(mucID)
00306         << " strip " << MucID::channel(mucID)
00307         << " Track Id " << (*iter_MucMcHit)->getTrackIndex() 
00308         << " pos x " << (*iter_MucMcHit)->getPositionX()
00309         << " pos y " << (*iter_MucMcHit)->getPositionY()
00310         << " pos z " << (*iter_MucMcHit)->getPositionZ()
00311         << endreq;
00312 
00313       McParticle *assocMcPart = 0;
00314       for (iter_mc = mcParticleCol->begin(); iter_mc != mcParticleCol->end(); iter_mc++) {
00315         if ( (*iter_mc)->trackIndex() == (*iter_MucMcHit)->getTrackIndex() ) {
00316           assocMcPart = *iter_mc;
00317           break;
00318         }
00319       }
00320       if (assocMcPart == 0) {
00321         log << MSG::WARNING << "No Corresponding Mc Particle found for this MucMcHit" << endreq;
00322       }
00323 
00324       MucMcHit *assocMucMcHit = *iter_MucMcHit;
00325       McPartToMucHitRel *relMcMuc = 0;
00326       relMcMuc = new McPartToMucHitRel( assocMcPart, assocMucMcHit );
00327       if (relMcMuc == 0) log << MSG::DEBUG << "relMcMuc not created " << endreq;
00328       else {
00329         bool addstat = assocMcMuc.addRelation( relMcMuc );
00330         if(!addstat) delete relMcMuc;
00331       }
00332     }
00333 
00334     log << MSG::DEBUG << " Fill McPartToMucHitTab, size " << assocMcMuc.size() << endreq;
00335     iter_mc = mcParticleCol->begin();
00336     for (;iter_mc != mcParticleCol->end(); iter_mc++) {
00337       log << MSG::DEBUG << " Track index " << (*iter_mc)->trackIndex() 
00338         << " particleId = " << (*iter_mc)->particleProperty()
00339         << endreq;   
00340 
00341       vector< McPartToMucHitRel* > vecMucMcHit = assocMcMuc.getRelByFirst(*iter_mc);
00342       vector< McPartToMucHitRel* >::iterator iter_MucMcHit = vecMucMcHit.begin();
00343       for (; iter_MucMcHit != vecMucMcHit.end(); iter_MucMcHit++) {
00344         mucID = (*iter_MucMcHit)->getSecond()->identify();
00345 
00346         log << MSG::DEBUG 
00347           //<< " McPartToMucHitTab " << iter_assocMcMuc
00348           << " MC Particle index " 
00349           << (*iter_MucMcHit)->getFirst()->trackIndex()
00350           << " contains " << " MucMcHit " 
00351           << " part "  << MucID::barrel_ec(mucID)
00352           << " seg "   << MucID::segment(mucID)
00353           << " gap "   << MucID::layer(mucID)
00354           << " strip " << MucID::channel(mucID)
00355           //<< " posX " << (*iter_MucMcHit)->getSecond()->getPositionX()
00356           //<< " posY " << (*iter_MucMcHit)->getSecond()->getPositionY()
00357           //<< " posZ " << (*iter_MucMcHit)->getSecond()->getPositionZ()
00358           << endreq;
00359       }
00360     }
00361 
00362     //assocMcPart = new McParticle(**iter_mc);
00363     //MucMcHit *assocMucMcHit = new MucMcHit(mucID, (*iter_MucMcHit)->getTrackIndex(), (*iter_MucMcHit)->getPositionX(),
00364     //                           (*iter_MucMcHit)->getPositionY(), (*iter_MucMcHit)->getPositionZ(),
00365     //                           (*iter_MucMcHit)->getPx(), (*iter_MucMcHit)->getPy(), (*iter_MucMcHit)->getPz() );
00366 
00367     // Retrieve McPartToMucHitTab 
00368     //SmartDataPtr<McPartToMucHitTab> aMcPartToMucHitTab(eventSvc(),"/Event/MC/McPartToMucHitTab");
00369     //if (!aMcPartToMucHitTab) {
00370     //log << MSG::ERROR << "Could not find McPartToMucHitTab" << endreq;
00371     //  return( StatusCode::FAILURE);
00372     //}
00373   }
00374 
00375   //
00376   // Read in Muc Digi;
00377   if (!m_MucRecHitContainer) {
00378     m_MucRecHitContainer = new MucRecHitContainer();
00379   }
00380   m_MucRecHitContainer->Clear();
00381   MucRecHitCol *aMucRecHitCol = new MucRecHitCol();
00382   m_MucRecHitContainer->SetMucRecHitCol(aMucRecHitCol);
00383 
00384   IDataManagerSvc *dataManSvc = dynamic_cast<IDataManagerSvc*> (eventSvc());
00385   DataObject *mucRecHitCol;
00386   eventSvc()->findObject("/Event/Recon/MucRecHitCol",mucRecHitCol);
00387   if(mucRecHitCol != NULL) {
00388     dataManSvc->clearSubTree("/Event/Recon/MucRecHitCol");
00389     eventSvc()->unregisterObject("/Event/Recon/MucRecHitCol");
00390   }
00391 
00392   sc = eventSvc()->registerObject("/Event/Recon/MucRecHitCol", aMucRecHitCol); 
00393   if (!sc) {
00394     log << MSG::ERROR << "/Event/Recon/MucRecHitCol not registerd!" << endreq;
00395     return( StatusCode::FAILURE);
00396   }
00397 
00398   log << MSG::INFO << "Add digis" << endreq; 
00399   MucDigiCol::iterator iter_Muc = mucDigiCol->begin();
00400   int mucDigiId = 0;
00401   for (;iter_Muc != mucDigiCol->end(); iter_Muc++, mucDigiId++) {
00402     mucID = (*iter_Muc)->identify();
00403     int part  = MucID::barrel_ec(mucID);
00404     int seg   = MucID::segment(mucID);
00405     int gap   = MucID::layer(mucID);
00406     int strip = MucID::channel(mucID);
00407     //m_MucRecHitContainer->AddHit(mucID);
00408     m_MucRecHitContainer->AddHit(part, seg, gap, strip);
00409 
00410     log << MSG::DEBUG <<  " digit" << mucDigiId << " : " 
00411       << " part "  << part
00412       << " seg "   << seg
00413       << " gap "   << gap
00414       << " strip " << strip
00415       << endreq;
00416   }
00417 
00418   //
00419   // Create track Container;
00420   RecMucTrackCol *aRecMucTrackCol = new RecMucTrackCol();
00421 
00422   //Register RecMucTrackCol to TDS
00423   DataObject *aReconEvent ;
00424   eventSvc()->findObject("/Event/Recon",aReconEvent);
00425   if(aReconEvent==NULL) {
00426     //then register Recon
00427     aReconEvent = new ReconEvent();
00428     StatusCode sc = eventSvc()->registerObject("/Event/Recon",aReconEvent);
00429     if(sc!=StatusCode::SUCCESS) {
00430       log << MSG::FATAL << "Could not register ReconEvent" <<endreq;
00431       return( StatusCode::FAILURE);
00432     }
00433   }
00434   StatusCode fr = eventSvc()->findObject("/Event/Recon",aReconEvent);
00435   if(fr!=StatusCode::SUCCESS) {
00436     log << MSG::WARNING << "Could not find register ReconEvent, will create it" <<endreq;          
00437     StatusCode sc = eventSvc()->registerObject("/Event/Recon",aReconEvent);
00438     if(sc!=StatusCode::SUCCESS) {
00439       log << MSG::FATAL << "Could not register ReconEvent" <<endreq;
00440       return( StatusCode::FAILURE);
00441     }
00442   }
00443 
00444   DataObject *mucTrackCol;
00445   eventSvc()->findObject("/Event/Recon/RecMucTrackCol",mucTrackCol);
00446   if(mucTrackCol != NULL) {
00447     dataManSvc->clearSubTree("/Event/Recon/RecMucTrackCol");
00448     eventSvc()->unregisterObject("/Event/Recon/RecMucTrackCol");
00449   }
00450 
00451   sc = eventSvc()->registerObject("/Event/Recon/RecMucTrackCol", aRecMucTrackCol);
00452   if(sc!=StatusCode::SUCCESS) {
00453     log << MSG::FATAL << "Could not register MUC track collection" << endreq;
00454     return( StatusCode::FAILURE);
00455   }
00456 
00457   // check RecMucTrackCol registered
00458   SmartDataPtr<RecMucTrackCol> findRecMucTrackCol(eventSvc(),"/Event/Recon/RecMucTrackCol");
00459   if (!findRecMucTrackCol) { 
00460     log << MSG::FATAL << "Could not find RecMucTrackCol" << endreq;
00461     return( StatusCode::FAILURE);
00462   }
00463   aRecMucTrackCol->clear();
00464 
00465   // m_ExtTrackSeedMode  1: ext from MC track, 2: use Ext track
00466 
00467   // Retrieve Ext track Col from TDS 
00468   //Check ExtTrackCol in TDS.
00469   SmartDataPtr<RecExtTrackCol> aExtTrackCol(eventSvc(),"/Event/Recon/RecExtTrackCol");
00470   if (!aExtTrackCol) {
00471     log << MSG::WARNING << "Can't find ExtTrackCol in TDS!" << endreq;
00472     //return( StatusCode::FAILURE);
00473   }
00474 
00475   SmartDataPtr<RecMdcTrackCol> aMdcTrackCol(eventSvc(),"/Event/Recon/RecMdcTrackCol");
00476   if (!aMdcTrackCol) {
00477     log << MSG::WARNING << "Can't find MdcTrackCol in TDS!" << endreq;
00478     //return( StatusCode::FAILURE);
00479   }
00480 
00481   //Retrieve Emc track col from TDS   2006.11.11 liangyt
00482   //Check EmcRecShowerCol in TDS.   
00483   SmartDataPtr<RecEmcShowerCol> aShowerCol(eventSvc(),"/Event/Recon/RecEmcShowerCol");
00484   if (!aShowerCol) { 
00485     log << MSG::WARNING << "Could not find RecEmcShowerCol" << endreq;
00486     //return( StatusCode::FAILURE);
00487   } 
00488 
00489   //    EmcRecShowerCol::iterator iShowerCol;
00490   //    for(iShowerCol=aShowerCol->begin();
00491   //        iShowerCol!=aShowerCol->end();
00492   //        iShowerCol++){
00493   //      cout<<"check EmcRecShowerCol registered:"<<endl;
00494   //          <<"shower position: "<<(*iShowerCol)->Position()<<endl;
00495   //    }
00496 
00497   int muctrackId = 0;
00498 
00499   log << MSG::INFO << "Add tracks info, ExtTrackSeedMode = " << m_ExtTrackSeedMode << ", EmcShowerSeedMode = " << m_EmcShowerSeedMode << endreq;
00500   //  if (m_ExtTrackSeedMode == 1 || !aExtTrackCol) {
00501   if (m_ExtTrackSeedMode == 1) 
00502   {  
00503     SmartDataPtr<McParticleCol> mcParticleCol(eventSvc(),"/Event/MC/McParticleCol");
00504     if (!mcParticleCol) {
00505       log << MSG::FATAL << "Could not find McParticle" << endreq;
00506       //return( StatusCode::FAILURE);
00507       return( StatusCode::SUCCESS);
00508     }
00509     McParticleCol::iterator iter_mc = mcParticleCol->begin();
00510 
00511     int trackIndex = -99; 
00512     for (int iTrack = 0;iter_mc != mcParticleCol->end(); iter_mc++, iTrack++) {
00513       if(!(*iter_mc)->primaryParticle()) continue;
00514       int pid = (*iter_mc)->particleProperty();
00515       int charge = 0;
00516       double mass = -1.0;
00517       if( pid >0 ) 
00518       { 
00519         if(m_particleTable->particle( pid )){
00520           charge = (int)m_particleTable->particle( pid )->charge();
00521           mass = m_particleTable->particle( pid )->mass();
00522         }
00523         else log << MSG::ERROR << "no this particle id in particle table, please check data" <<endreq;
00524       } else if ( pid <0 ) 
00525       {
00526         if(m_particleTable->particle( -pid )){
00527           charge = (int)m_particleTable->particle( -pid )->charge();
00528           charge *= -1;
00529           mass = m_particleTable->particle( -pid )->mass();
00530         }
00531         else log << MSG::ERROR << "no this particle id in particle table, please check data" <<endreq;
00532       } else {
00533         log << MSG::WARNING << "wrong particle id, please check data" <<endreq;
00534       }
00535 
00536       if(!charge) {
00537         log << MSG::WARNING
00538           << " neutral particle charge = 0!!! ...just skip it !" 
00539           << endreq;
00540         continue;
00541       }
00542       
00543       trackIndex = (*iter_mc)->trackIndex();
00544       log << MSG::DEBUG << "iTrack " << iTrack << " index " << trackIndex
00545                 << " particleId = " << (*iter_mc)->particleProperty() 
00546                 << endreq;   
00547       
00548       RecMucTrack *aTrack = new RecMucTrack();
00549       aTrack->setTrackId(trackIndex);
00550       aTrack->setId(muctrackId);
00551       
00552       HepLorentzVector initialMomentum = (*iter_mc)->initialFourMomentum();
00553       HepLorentzVector initialPos = (*iter_mc)->initialPosition();
00554       float theta0 = initialMomentum.theta();
00555       float phi0   = initialMomentum.phi();
00556       //float dirX = sin(theta0) * cos(phi0);
00557       //float dirY = sin(theta0) * sin(phi0);
00558       //float dirZ = cos(theta0);
00559       float x0 = initialPos.x();
00560       float y0 = initialPos.y();
00561       float z0 = initialPos.z();
00562       
00563       Hep3Vector startPos(x0, y0, z0);
00564       Hep3Vector startP(initialMomentum.px(), initialMomentum.py(), initialMomentum.pz());
00565       log << MSG::DEBUG << "startP " << startP <<"  startPos "<<startPos<< endreq;
00566       Hep3Vector endPos(0, 0, 0), endP;
00567       
00568       aTrack->GetMdcExtTrack(startPos, startP, charge, endPos, endP);
00569       aTrack->SetMdcPos( x0, y0, z0);
00570       aTrack->SetMdcMomentum( startP.x(), startP.y(), startP.z() );
00571       aTrack->SetExtTrackID(trackIndex);
00572       aTrack->SetExtMucPos( endPos.x(), endPos.y(), endPos.z() );
00573       aTrack->SetExtMucMomentum( endP.x(), endP.y(), endP.z() );
00574       aTrack->SetMucPos( endPos.x(), endPos.y(), endPos.z() );
00575       aTrack->SetMucMomentum( endP.x(), endP.y(), endP.z() );
00576       aTrack->SetCurrentPos( endPos.x(), endPos.y(), endPos.z());
00577       aTrack->SetCurrentDir( endP.x(), endP.y(), endP.z() );
00578       //aTrack->SetMucVertexPos( x0, y0, z0);
00579       //aTrack->SetMucVertexDir( dirX, dirY, dirZ );
00580       //aTrack->SetCurrentPos( x0, y0, z0);
00581       //aTrack->SetCurrentDir( dirX, dirY, dirZ );
00582       
00583       //log << MSG::DEBUG 
00584       //<< "  pdg " << (*iter_mc)->particleProperty()
00585       //<< "  mucPos " << aTrack->GetMucVertexPos() 
00586       //<< "  mucDir " << aTrack->GetMucVertexDir() 
00587       //<< endreq; 
00588       log << MSG::INFO << "Add track, trkIndex: " << trackIndex << "\tmucTrkId: " << muctrackId << endreq; 
00589       aRecMucTrackCol->add(aTrack);
00590       muctrackId ++ ;
00591     }
00592   }
00593   else if (m_ExtTrackSeedMode == 2) 
00594   {
00595     if (!aExtTrackCol || !aMdcTrackCol) {
00596       log << MSG::WARNING << "Can't find ExtTrackCol or MdcTrackCol in TDS!" << endreq;
00597       return StatusCode::SUCCESS;
00598     }
00599 
00600     int trackIndex = -99; 
00601     //log << MSG::INFO << "MdcTrackCol:\t " << aMdcTrackCol->size() << "\tExtTrackCol:\t " << aExtTrackCol->size() << endreq; 
00602     RecExtTrackCol::iterator iter_ExtTrack = aExtTrackCol->begin();
00603     RecMdcTrackCol::iterator iter_MdcTrack = aMdcTrackCol->begin();
00604     int iExtTrack = 0;
00605     for (;iter_ExtTrack != aExtTrackCol->end() && iter_MdcTrack != aMdcTrackCol->end(); iter_ExtTrack++, iter_MdcTrack++, iExtTrack++) 
00606     { 
00607       trackIndex = (*iter_ExtTrack)->GetTrackId();
00608     log << MSG::DEBUG << "iExtTrack " << iExtTrack << " index " << trackIndex << " MucPos "      
00609         << iExtTrack << (*iter_ExtTrack)->mucPosition().x() << " " 
00610         << (*iter_ExtTrack)->mucPosition().y() << " " 
00611         << (*iter_ExtTrack)->mucPosition().z() << " r " 
00612         << (*iter_ExtTrack)->mucPosition().r() << endreq;
00613 
00614       if((*iter_ExtTrack)->mucPosition().x() == -99 &&
00615           (*iter_ExtTrack)->mucPosition().y() == -99 &&
00616           (*iter_ExtTrack)->mucPosition().z() == -99)
00617       {
00618         log << MSG::INFO <<"Bad ExtTrack, trackIndex: " << trackIndex << ", skip" << endreq;
00619         continue;
00620       }
00621 
00622       RecMucTrack *aTrack = new RecMucTrack();
00623       aTrack->setTrackId(trackIndex);
00624       aTrack->setId(muctrackId);
00625 
00626       aTrack->setType(0); //0 for use seed from mdc ext;
00627       
00628       aTrack->SetExtTrack(*iter_ExtTrack);
00629 
00630       //Hep3Vector mdcPos      = (*iter_ExtTrack)->GetMdcPosition();
00631       Hep3Vector mdcMomentum = (*iter_MdcTrack)->p3();
00632       Hep3Vector mucPos      = (*iter_ExtTrack)->mucPosition();
00633       Hep3Vector mucMomentum = (*iter_ExtTrack)->mucMomentum();
00634     
00635       //aTrack->SetMdcPos( mdcPos.x(), mdcPos.y(), mdcPos.z() );
00636       aTrack->SetMdcMomentum( mdcMomentum.x(), mdcMomentum.y(), mdcMomentum.z() );
00637       aTrack->SetExtMucPos( mucPos.x(), mucPos.y(), mucPos.z() );
00638       //cout << "ExtTrack MucPos " << aTrack->GetExtMucPos() << endl;
00639       aTrack->SetExtMucMomentum( mucMomentum.x(), mucMomentum.y(), mucMomentum.z() );
00640       aTrack->SetMucPos( mucPos.x(), mucPos.y(), mucPos.z() );
00641       aTrack->SetMucMomentum( mucMomentum.x(), mucMomentum.y(), mucMomentum.z() );
00642       aTrack->SetCurrentPos( mucPos.x(), mucPos.y(), mucPos.z());
00643       aTrack->SetCurrentDir( mucMomentum.x(), mucMomentum.y(), mucMomentum.z() );
00644       aTrack->SetRecMode(0);  //mdc ext mode
00645       
00646       log << MSG::INFO << "Add track, trkIndex: " << trackIndex << "\tmucTrkId: " << muctrackId << endreq; 
00647       aRecMucTrackCol->add(aTrack);
00648       muctrackId ++ ;
00649     }
00650   }
00651   else if(m_ExtTrackSeedMode == 3) 
00652   { //cosmic ray
00653 
00654     if (!aMdcTrackCol) {
00655       log << MSG::WARNING << "Can't find MdcTrackCol in TDS!" << endreq;
00656       return StatusCode::SUCCESS;
00657       //return( StatusCode::FAILURE);
00658     }
00659     log<< MSG::INFO << "Mdc track size: "<<aMdcTrackCol->size()<<endreq;
00660 
00661     int trackIndex = -99; 
00662     for(RecMdcTrackCol::iterator iter_mdc1 = aMdcTrackCol->begin(); iter_mdc1 != aMdcTrackCol->end(); iter_mdc1++){
00663       //if((*iter_mdc1)->getFi0() > 0.5*pi && (*iter_mdc1)->getFi0() < 1.5*pi){
00664       //  cout<<"this is down track"<<endl;
00665       //}
00666     
00667       trackIndex = (*iter_mdc1)->trackId();
00668       HepVector helix = (*iter_mdc1)->helix();
00669       //cout<<"helix: "<<helix<<endl;
00670 
00671       float x0 = helix[0] * cos(helix[1]);
00672       float y0 = helix[0] * sin(helix[1]);
00673       float z0 = helix[3];
00674 
00675       //float dx = 1/(sqrt(1+helix[4]*helix[4])) * (-1* sin(helix[1]));
00676       //float dy = 1/(sqrt(1+helix[4]*helix[4])) * cos(helix[1]);
00677       //float dz = 1/(sqrt(1+helix[4]*helix[4])) * helix[4];
00678   
00679       float dx =  -1* sin(helix[1]);
00680       float dy =  cos(helix[1]);
00681       float dz =  helix[4];
00682 
00683       //cout<<"pos: "<<x0<<"  "<<y0<<"  "<<z0<<"  dir: "<<dx<<"  "<<dy<<"  "<<dz<<endl;
00684 
00685       RecMucTrack *aTrack = new RecMucTrack();
00686       aTrack->setTrackId(trackIndex);
00687       aTrack->setId(muctrackId);
00688       muctrackId ++ ;
00689 
00690       aTrack->setType(3); //3 for use seed from mdc without magnet field;
00691 
00692     
00693       Hep3Vector mucPos(x0,y0,z0);
00694       Hep3Vector mucMomentum(dx,dy,dz);
00695     
00696       aTrack->SetExtMucPos(x0,y0,z0);
00697       aTrack->SetExtMucMomentum(dx,dy,dz);
00698     
00699       aTrack->SetMucPos(x0,y0,z0);
00700       aTrack->SetMucMomentum(dx,dy,dz);
00701       aTrack->SetCurrentPos(x0,y0,z0);
00702       aTrack->SetCurrentDir(dx,dy,dz);
00703       aTrack->SetRecMode(0);  //mdc ext mode
00704       
00705       log << MSG::INFO << "Add track, trkIndex: " << trackIndex << "\tmucTrkId: " << muctrackId << endreq; 
00706       aRecMucTrackCol->add(aTrack);
00707       muctrackId ++ ;
00708     }
00709   }
00710   else {log << MSG::INFO <<"ExtTrackSeedMode error!"<<endreq; }
00711 
00712   //Rec from Emc extend track
00713   if(m_EmcShowerSeedMode == 1)
00714   {
00715     int trackIndex = 999;
00716     RecEmcShowerCol::iterator iShowerCol;
00717     for(iShowerCol=aShowerCol->begin();  iShowerCol!=aShowerCol->end(); iShowerCol++)
00718     {
00719       //       cout<<"shower id: "<<(*iShowerCol)->ShowerId()<<" "
00720       //          <<"shower energy: "<<(*iShowerCol)->Energy()<<" "
00721       //          <<"shower position: "<<(*iShowerCol)->Position()<<endl;
00722 
00723       RecMucTrack *aTrack = new RecMucTrack();
00724       aTrack->setTrackId(trackIndex);
00725       aTrack->setId(muctrackId);
00726       aTrack->setType(1); //1 for use seed from emc hits;
00727 
00728       Hep3Vector mucPos      = (*iShowerCol)->position();
00729       Hep3Vector mucMomentum = (*iShowerCol)->position();
00730 
00731       aTrack->SetExtMucPos( mucPos.x(), mucPos.y(), mucPos.z() );
00732       //cout << "EmcTrack MucPos " << aTrack->GetExtMucPos() << endl;
00733       aTrack->SetExtMucMomentum( mucMomentum.x(), mucMomentum.y(), mucMomentum.z() );
00734       aTrack->SetMucPos( mucPos.x(), mucPos.y(), mucPos.z() );
00735       aTrack->SetMucMomentum( mucMomentum.x(), mucMomentum.y(), mucMomentum.z() );
00736       aTrack->SetCurrentPos( mucPos.x(), mucPos.y(), mucPos.z());
00737       aTrack->SetCurrentDir( mucMomentum.x(), mucMomentum.y(), mucMomentum.z() );
00738       aTrack->SetRecMode(1);  //emc ext mode
00739       //cout<<" MucRecTrkExt 111 trackIndex :"<<aTrack->GetIndex()<<endl;
00740       
00741       log << MSG::INFO << "Add track, trkIndex: " << trackIndex << "\tmucTrkId: " << muctrackId << endreq; 
00742       aRecMucTrackCol->add(aTrack);
00743       muctrackId ++ ;
00744 
00745       m_emcrec = 1;
00746     }
00747   }
00748   log << MSG::DEBUG << " track container filled " << endreq;
00749 
00750   // All input filled, begin track finding;
00751   log << MSG::INFO << "Start tracking..." << endreq;
00752   int runVerbose = 1;
00753   RecMucTrack *aTrack = 0;
00754   for(int iTrack = 0; iTrack < (int)aRecMucTrackCol->size(); iTrack++) 
00755   {
00756     log << MSG::DEBUG << "iTrack " << iTrack << endreq;
00757     aTrack = (*aRecMucTrackCol)[iTrack]; 
00758     //cout<<"in MucRecTrkExt trackIndex :"<<aTrack->GetIndex()<<endl;
00759     
00760     Hep3Vector currentPos = aTrack->GetCurrentPos();
00761     Hep3Vector currentDir = aTrack->GetCurrentDir();
00762     if(currentPos.mag() < kMinor) {
00763       log << MSG::WARNING << "No MUC intersection in track " << iTrack << endreq;
00764       continue;
00765     }
00766     //log << MSG::INFO << " pos " << currentPos << " dir " << currentDir << endreq;
00767     
00768     int firstHitFound[2] = { 0,  0}; // Has the fist position in this orient determined? if so, could CorrectDirection()
00769     int firstHitGap[2]   = {-1, -1}; // When the fist position in this orient determined, the gap it is on
00770     for(int partSeq = 0; partSeq < (int)MucID::getPartNum(); partSeq++) 
00771     {
00772       int iPart = kPartSeq[partSeq];
00773       for(int iGap = 0; iGap < (int)MucID::getGapNum(iPart); iGap++) 
00774       {
00775         //log << MSG::INFO << iPart << iGap << endreq;
00776         int seg = -1;
00777         int orient = MucGeoGeneral::Instance()->GetGap(iPart, 0, iGap)->Orient();;
00778 
00779         float xInsct, yInsct, zInsct;
00780         aTrack->Project(iPart, iGap, xInsct, yInsct, zInsct, seg);
00781         if(m_MsOutput) cout << "part " << iPart << " gap " << iGap << " x " << xInsct << " y " << yInsct << " z " << zInsct << " seg " << seg << endl; 
00782         
00783         if(seg == -1) continue;
00784 
00785         aTrack->SetCurrentInsct(xInsct, yInsct, zInsct);
00786 
00787         for(int iDeltaSeg = 0; iDeltaSeg < kNSegSearch; iDeltaSeg++) 
00788         {
00789           int iSeg = seg + kDeltaSeg[iDeltaSeg];  // also find on neighbor seg;
00790           if(iSeg < 0) iSeg += MucID::getSegNum(iPart);
00791           if(iSeg > (int)MucID::getSegNum(iPart) - 1) iSeg -= MucID::getSegNum(iPart);
00792 
00793           //log << MSG::DEBUG << iPart << iSeg << iGap 
00794           //  << "NHits " << m_MucRecHitContainer->GetGapHitCount(iPart, seg, iGap) << endreq;
00795 
00796           //----------now change window size(depond on mom)
00797           //float window = kWindowSize[iPart][iGap];
00798           Hep3Vector mom_mdc = aTrack->getMdcMomentum();
00799           float window = getWindowSize(mom_mdc, iPart, iSeg, iGap);
00800 
00801           if (firstHitFound[orient] != 1) window *= kFirstHitWindowRatio;   // if no hits have been found on this orient, expand the window.
00802           for (int iHit = 0; iHit < m_MucRecHitContainer->GetGapHitCount(iPart, iSeg, iGap); iHit++) 
00803           {
00804             log << MSG::DEBUG << "iSeg " << iSeg << " iHit " << iHit << endreq;
00805             MucRecHit* pHit = m_MucRecHitContainer->GetHit(iPart, iSeg, iGap, iHit);
00806             //cout<< "strip " << pHit->Strip() << " center " << pHit->GetCenterPos() << endl;
00807 
00808             if (!pHit) {
00809               log << MSG::WARNING << "MucRecTrkExt: null pointer to pHit" << endreq;
00810               continue;
00811             }
00812             else 
00813             {
00814               //cout<<"in MucRecTrkExt: pHit exist : "<<iPart<<" "<<iSeg<<" "<<iGap<<" "<<iHit<<" mode:"<<pHit->GetHitMode()<<endl;
00815               // Get the displacement of hit pHit to intersection
00816               //float dX  = aTrack->GetHitDistance(pHit);
00817               float dX  = aTrack->GetHitDistance2(pHit);  //not abs value now!
00818               log << MSG::DEBUG << "distance = " << setw(8) << dX << " size " << setw(4) << window << endreq;
00819 
00820               if(m_NtOutput >= 2){ //too many info
00821                 m_distance = dX;
00822                 m_part = iPart; m_seg = iSeg; m_gap = iGap; m_strip = pHit->Strip();
00823                 m_diff = -99; m_tuple->write();
00824               }
00825 
00826               if (fabs(dX) < window) 
00827               {
00828                 // Attach the hit if it exists
00829                 // cout << "meet window " << pHit->GetSoftID() << endl;
00830                 //****************if this if emc track, we abondon used hit in mdc*******************
00831                 //if(m_emcrec == 1 )
00832                 if(aTrack->GetRecMode() == 0) {
00833                   pHit->SetHitMode(1); //mdc ext
00834                   aTrack->AttachHit(pHit);
00835                   //cout<<"in MucRecTrkExt: trackmode==0 so mdc ext "<<iPart<<" "<<iSeg<<" "<<iGap<<" "<<iHit<<endl;
00836                 }
00837                 if(aTrack->GetRecMode() == 1) {
00838                   //cout<<"in MucRecTrkExt: HitMode = "<<pHit->GetHitMode()<<" part:  "<<iPart<<" "<<iSeg<<" "<<iGap<<" "<<iHit<<endl;
00839                   if(pHit->GetHitMode()!=1) {
00840                     aTrack->AttachHit(pHit);  //this hit has not been used by mdc ext
00841                     pHit->SetHitMode(2); //emc ext
00842                   }
00843                 }
00844 
00845                 //push back distance between ext and hits
00846                 aTrack->pushExtDistHits(dX);   //2009-03-12
00847 
00848                 if (firstHitGap[orient] == -1) firstHitGap[orient] = iGap;
00849                 firstHitFound[orient] = 1;
00850                 //cout << "You could correct directon in orient " << orient << endl;
00851 
00852                 log << MSG::DEBUG << " part " << iPart << " seg " << iSeg << " gap " << iGap 
00853                     << " strip " << setw(2) << pHit->Strip() << " attatched" << endreq;
00854                 log << MSG::DEBUG << "current total hits " << aTrack->GetTotalHits() << endreq;
00855               }
00856               else {
00857                 m_NHitsLostInGap[iGap]++;
00858                 //log << MSG::DEBUG cout << " part " << iPart << " seg " << iSeg << " gap " << iGap 
00859                 //       << " strip " << setw(2) << pHit->GetSoftID().GetStrip() 
00860                 //     << " not attached !" << endreq;
00861               }
00862             } // end pHit else 
00863           } // end for iHit
00864           aTrack->CalculateInsct(iPart, iSeg, iGap);
00865         } // end for iDeltaSeg
00866         
00867         // When correct dir in the orient is permitted and this gap is not the gap first hit locates.
00868         if(m_ExtTrackSeedMode != 3 && !m_Blind) { //for cosmic ray, we do not correct dir and pos...
00869           if (firstHitFound[orient] && firstHitGap[orient] != iGap) aTrack->CorrectDir();
00870           aTrack->CorrectPos();
00871         }
00872         //cout << "Current Intersection " << aTrack->GetCurrentInsct() << endl;
00873         //cout << "Current Direction    " << aTrack->GetCurrentDir() << endl;
00874         aTrack->AttachInsct(aTrack->GetCurrentInsct());
00875       } // end for iGap
00876     } // end for iSeg
00877     aTrack->LineFit(m_fittingMethod);
00878     aTrack->ComputeTrackInfo(m_fittingMethod);
00879     log << MSG::INFO <<"Fit track done! trackIndex: " << aTrack->trackId() << ", mucId: " << aTrack->id() << ", RecMode: " << aTrack->GetRecMode() << endreq;
00880  
00881     if(m_NtOutput>=1)
00882     {
00883       m_depth = aTrack->depth();
00884       m_Distance = aTrack->distance();
00885       m_MaxHits = aTrack->maxHitsInLayer();
00886       m_Chi2  = aTrack->chi2();
00887       m_Dist_x = aTrack->GetExtMucDistance().x();
00888       m_Dist_y = aTrack->GetExtMucDistance().y();
00889       m_Dist_z = aTrack->GetExtMucDistance().z();
00890       m_posx_ext = aTrack->GetMucStripPos().x();
00891       m_posy_ext = aTrack->GetMucStripPos().y();
00892       m_posz_ext = aTrack->GetMucStripPos().z();
00893 
00894       m_emctrack = m_emcrec;
00895     }
00896 
00897     //test distance of line/quad fitting
00898     if(m_NtOutput>=2)
00899     {
00900             vector<MucRecHit*> attachedHits = aTrack->GetHits();
00901             vector<MucRecHit*> expectedHits = aTrack->GetExpectedHits();
00902             vector<float>      distanceHits = aTrack->getDistHits();
00903             vector<float>      distanceHits_quad =  aTrack->getQuadDistHits();
00904             vector<float>      distanceHits_ext =  aTrack->getExtDistHits();
00905 
00906             for(int i=0; i< expectedHits.size(); i++)
00907             {
00908                     MucRecHit *ihit = expectedHits[i];
00909                     //cout<<"expected Hits: "<<ihit->Part()<<"  "<<ihit->Seg()<<"  "<<ihit->Gap()<<"  "<<ihit->Strip()<<endl;
00910             }
00911 
00912             for(int j=0; j< attachedHits.size(); j++)
00913       {
00914                     MucRecHit *jhit = attachedHits[j];
00915                     //                if(attachedHits.size()!=distanceHits.size())cout<<"attached hits size no same as disthits!!!"<<endl;
00916                     if(attachedHits.size()==distanceHits.size())
00917                     {    // same gap, cale distance
00918                             m_part = jhit->Part();   m_seg = jhit->Seg();  m_gap = jhit->Gap();
00919                             m_strip = jhit->Strip();
00920                             m_distance = distanceHits[j];
00921                             m_Distance = distanceHits_quad[j];
00922                             m_diff = -9999;
00923                             //cout<<"distance = "<<m_dist<<endl;
00924                             //if(distanceHits.size() == distanceHits_ext.size())
00925                             //cout<<"ext dist: "<<distanceHits_ext[j]<<endl;
00926                             m_tuple->write();
00927                     }
00928             }
00929     } // end if output
00930 
00931     int nHitsAttached = aTrack->GetTotalHits();
00932     //m_NHitsAttachedTotal += nHitsAttached;
00933     //if(mucDigiCol->size() - recHitNum != 0) 
00934     log << MSG::DEBUG << "track Index " << aTrack->trackId()
00935               << setw(2) << mucDigiCol->size() - nHitsAttached << " of "
00936               << setw(2) << mucDigiCol->size() << " lost " << endreq;
00937     //m_NHitsLost[int(mucDigiCol->size() - nHitsAttached)]++;
00938   } // end for iTrack
00939 
00940   
00941   // if (aRecMucTrackCol->size() == 0) m_NHitsLost[0]++;
00942   //m_NHitsTotal         += mucDigiCol->size();
00943     
00944   //****************** look up unrec hits to do recon with MUC info ***************??
00945 
00946   if(m_MucHitSeedMode == 1)
00947   {
00948     MucRecHit *pHit = 0 ,*pHit0 = 0, *pHit1 = 0;
00949     int count, orient;
00950 
00951     for (int iPart = 0; iPart < (int)MucID::getPartNum(); iPart++) 
00952     {
00953       for (int iSeg = 0; iSeg < (int)MucID::getSegNum(iPart); iSeg++) 
00954       {
00955         bool hit0 = false, hit1 = false;  int firstgap0 = -1, firstgap1 = -1; int nStrip0 = 0, nStrip1 = 0;
00956         Hep3Vector posHit0, posHit1;
00957         for(int iGap = 0; iGap < (int)MucID::getGapNum(iPart); iGap++) 
00958         {
00959           count = m_MucRecHitContainer->GetGapHitCount(iPart, iSeg, iGap);
00960           for (int iHit0 = 0; iHit0 < count; iHit0++) 
00961           {
00962             //cout << "iHit0 = " << iHit0 << endl;
00963             pHit = m_MucRecHitContainer->GetHit(iPart, iSeg, iGap, iHit0);
00964             if (!pHit) {
00965               log << MSG::FATAL << "MucRecRoadFinder-E10 " << " part " << iPart
00966                 << " seg " << iSeg << " gap "  << iGap << " null pointer"
00967                 << endl;
00968             }
00969             else 
00970             {
00971               //cout<< "pHit Mode is : " << pHit->GetHitMode() << endl;
00972               if(pHit->GetHitMode() == -1)
00973               {
00974                 orient = MucGeoGeneral::Instance()->GetGap(iPart, 0, iGap)->Orient();
00975                 if(orient == 1 && hit0 == false){
00976                   hit0 = true;
00977                   firstgap0 = iGap;
00978                 }
00979                 if(iGap == firstgap0){
00980                   posHit0 += pHit->GetCenterPos();
00981                   nStrip0++;
00982                 }
00983 
00984                 if(orient == 0 && hit1 == false){
00985                   hit1 = true;
00986                   firstgap1 = iGap;
00987                 }
00988                 if(iGap == firstgap1){
00989                   posHit1 += pHit->GetCenterPos();
00990                   nStrip1++;
00991                 }
00992                 //            if(orient == 1 && hit0 == false){   //orient = 1
00993                 //              hit0 = true;
00994                 //              pHit0 = pHit;  //pHit0 is the first hit of first gap in this segment with orient 1
00995                 //            }
00996                 //            if(orient == 0 && hit1 == false){
00997                 //              hit1 = true;
00998                 //              pHit1 = pHit;  //pHit0 is the first hit of first gap in this segment with orient 0
00999                 //            }
01000               }
01001             }// pHit0 exist;
01002           } // iHit0
01003         }  //iGap
01004 
01005         //if hit0 hit1 exist, make a seed
01006         if(hit0 && hit1)
01007         {
01008           posHit0 /= nStrip0;
01009           posHit1 /= nStrip1;
01010           //cout<< "pHit0 position is " << posHit0  <<" "<<nStrip0<<endl;
01011           //cout<< "pHit1 position is " << posHit1  <<" "<<nStrip1<<endl;
01012 
01013           int trackIndex = 999;
01014           RecMucTrack *aTrack = new RecMucTrack();
01015           aTrack->setTrackId(trackIndex);
01016           aTrack->setId(muctrackId);
01017           muctrackId ++ ;
01018 
01019           aTrack->setType(2); //2 for use seed from Muc Information
01020 
01021           Hep3Vector mucPos, mucMomentum;    
01022           if(iPart==1) {
01023             mucMomentum.set(posHit0.x(),posHit0.y(),posHit1.z());
01024             mucPos = mucMomentum * 0.9;  //shorten mucPos, otherwise maybe no intersection point found for the track and this gap
01025           }
01026           else{
01027             mucMomentum.set(posHit0.x(),posHit1.y(),posHit0.z()*0.5+posHit1.z()*0.5);
01028             mucPos = mucMomentum * 0.9;  //shorten mucPos, otherwise maybe no intersection point found for the track and this gap
01029           }
01030 
01031           aTrack->SetExtMucPos( mucPos.x(), mucPos.y(), mucPos.z() );
01032           //cout << "EmcTrack MucPos " << aTrack->GetExtMucPos() << endl;
01033           aTrack->SetExtMucMomentum( mucMomentum.x(), mucMomentum.y(), mucMomentum.z() );
01034           aTrack->SetMucPos( mucPos.x(), mucPos.y(), mucPos.z() );
01035           aTrack->SetMucMomentum( mucMomentum.x(), mucMomentum.y(), mucMomentum.z() );
01036           aTrack->SetCurrentPos( mucPos.x(), mucPos.y(), mucPos.z());
01037           aTrack->SetCurrentDir( mucMomentum.x(), mucMomentum.y(), mucMomentum.z() );
01038           //cout<<" MucRecTrkExt 111 trackIndex :"<<aTrack->GetIndex()<<endl;
01039           aTrack->SetRecMode(2); 
01040           aRecMucTrackCol->add(aTrack);
01041 
01042           TrackFinding(aTrack);
01043         }
01044 
01045         if(m_NtOutput>=1){
01046           m_depth = aTrack->depth();
01047           m_Distance = aTrack->distance();
01048           m_MaxHits = aTrack->maxHitsInLayer();
01049           m_Chi2  = aTrack->chi2();
01050           m_Dist_x = aTrack->GetExtMucDistance().x();
01051           m_Dist_y = aTrack->GetExtMucDistance().y();
01052           m_Dist_z = aTrack->GetExtMucDistance().z();
01053           m_posx_ext = aTrack->GetMucStripPos().x();
01054           m_posy_ext = aTrack->GetMucStripPos().y();
01055           m_posz_ext = aTrack->GetMucStripPos().z();
01056 
01057           m_emctrack = 2;
01058         }
01059 
01060       }  //iSeg
01061     }  //iPart
01062   } // end if SeedMode == 1
01063 
01064   //************** end of look up unrec hits to do recon with MUC info ************??
01065   int nTracksTotal = 0;//mcParticleCol->size();
01066   int nTracksFound = 0;
01067   int nTracksLost  = 0;
01068   int nTracksLostByExt = 0;
01069   int nTracksMisFound = 0;
01070   
01071   int nDigisTotal   = 0;//mucDigiCol->size();
01072   int nHitsTotal    = 0;//assocMcMuc.size();
01073   int nHitsFound    = 0;
01074   int nHitsLost     = 0;
01075   int nHitsMisFound = 0;
01076   
01077 /*
01078   if (m_CompareWithMcHit) {
01079     // 
01080     // Compare RecMucTracks with McPartToMucHitTab;
01081 
01082     log << MSG::DEBUG << " *******************************" << endreq;
01083     log << MSG::DEBUG << " Compare Mc Info with rec tracks" << endreq;
01084     log << MSG::DEBUG << " McParticleCol     size " << mcParticleCol->size() << endreq;
01085     log << MSG::DEBUG << " McPartToMcuHitTab size " << assocMcMuc.size() << endreq;
01086 
01087     iter_mc = mcParticleCol->begin();
01088     for (;iter_mc != mcParticleCol->end(); iter_mc++) {
01089       log << MSG::DEBUG << " McParticle index " << (*iter_mc)->getTrackIndex() 
01090         << " particleId = " << (*iter_mc)->particleProperty()
01091         << endreq;   
01092 
01093       bool foundRecoTrack = false;
01094       log << MSG::DEBUG << " Reconed RecMucTrackCol size " << aRecMucTrackCol->size() << endreq;
01095       aTrack = 0;
01096       for(int iTrack = 0; iTrack < (int)aRecMucTrackCol->size(); iTrack++) {
01097         log << MSG::DEBUG << "iTrack " << iTrack << endreq;
01098         aTrack = (*aRecMucTrackCol)[iTrack]; 
01099         if (aTrack->GetExtTrackID() == (*iter_mc)->getTrackIndex()) {
01100           log << MSG::DEBUG << " found iTrack " << iTrack 
01101             << " index " << aTrack->GetExtTrackID() 
01102             << " equals to " << " McParticle index " << (*iter_mc)->getTrackIndex()
01103             << endreq;
01104           foundRecoTrack = true;
01105           break;
01106         }
01107       }
01108       if (foundRecoTrack == false) {
01109         nTracksLost++;
01110         m_TrackLostByMdc.push_back(eventHeader->eventNumber());
01111         log << MSG::DEBUG << " this McParticle is not found in RecMucTrackCol,"
01112           << "->not intialized from ExtTrack-> not constructed in Mdc" << endreq;
01113       }
01114       else {
01115         nTracksFound++;
01116       }
01117 
01118       int nHitsFoundInTrack   = 0; 
01119       int nHitsLostInTrack    = 0; 
01120 
01121       vector< McPartToMucHitRel* > vecMucMcHit = assocMcMuc.getRelByFirst(*iter_mc);
01122       log << MSG::DEBUG << " McParticle " << (*iter_mc)->getTrackIndex() 
01123         << " contains " << vecMucMcHit.size() << " MucMcHit " << endreq;
01124       if (!aTrack) {
01125         nHitsLostInTrack = vecMucMcHit.size();
01126         m_NHitsLostByMdcTotal += nHitsLostInTrack;
01127         log << MSG::DEBUG << " This Mc particle has no corresponding Reco Track, "
01128           << " all of its MucMcHits lost" << endreq;  
01129       }
01130       else if ((aTrack->GetExtMucPos()).mag() < 0.1) {
01131         nHitsLostInTrack = vecMucMcHit.size();
01132         m_NHitsLostByExtTotal += nHitsLostInTrack;
01133         nTracksLostByExt++;
01134         m_TrackLostByExt.push_back(eventHeader->eventNumber());
01135         log << MSG::DEBUG << " This Mc particle 's Ext track's intersection with Muc is zero," 
01136           << " all of its MucMcHits lost" << endreq;
01137       }
01138       else {
01139         vector< McPartToMucHitRel* >::iterator iter_MucMcHit = vecMucMcHit.begin();
01140         for (; iter_MucMcHit != vecMucMcHit.end(); iter_MucMcHit++) {
01141           mucID = (*iter_MucMcHit)->getSecond()->identify();
01142           int part  = MucID::barrel_ec(mucID);
01143           int seg   = MucID::segment(mucID);
01144           int gap   = MucID::layer(mucID);
01145           int strip = MucID::channel(mucID);
01146 
01147           log << MSG::DEBUG 
01148             //<< " McPartToMucHitTab " << iter_assocMcMuc
01149             << " McParticle index " 
01150             << (*iter_MucMcHit)->getFirst()->getTrackIndex()
01151             << " contains " << " MucMcHit " 
01152             << " part "  << part
01153             << " seg "   << seg
01154             << " gap "   << gap
01155             << " strip " << strip
01156             //<< " posX " << (*iter_MucMcHit)->getSecond()->getPositionX()
01157             // << " posY " << (*iter_MucMcHit)->getSecond()->getPositionY()
01158             //<< " posZ " << (*iter_MucMcHit)->getSecond()->getPositionZ()
01159             << endreq;
01160 
01161           if (aTrack->HasHit(part, seg, gap, strip)) {
01162             nHitsFoundInTrack++;
01163             log << MSG::DEBUG << " This MucMchit found on this Reco track" << endreq;
01164           }
01165           else { 
01166             nHitsLostInTrack++;
01167             log << MSG::DEBUG << " This MucMchit not found on this Reco track, it is lost " << endreq;
01168           }
01169         }
01170       }
01171 
01172       nHitsFound += nHitsFoundInTrack;
01173       nHitsLost  += nHitsLostInTrack;
01174 
01175       m_NHitsLost[nHitsLost]++;
01176       if (nHitsLost > 0) {
01177         m_TrackLostHit.push_back(eventHeader->eventNumber());
01178         m_TrackLostHitNumber.push_back(nHitsLost);
01179       }
01180 
01181       log << MSG::DEBUG << " In " << vecMucMcHit.size() << " MucMcHit : "
01182         << nHitsFoundInTrack << " found in Reco track, " 
01183         << nHitsLostInTrack << " lost " << endreq;
01184     }
01185 
01186     // Reversely, Compare McPartToMucHitTab with RecMucTracks;
01187     log << MSG::DEBUG << " *******************************" << endreq;
01188     log << MSG::DEBUG << " Compare rec tracks with Mc Info" << endreq;
01189     log << MSG::DEBUG << " Reconed RecMucTrackCol size " << aRecMucTrackCol->size() << endreq;
01190     for(int iTrack = 0; iTrack < (int)aRecMucTrackCol->size(); iTrack++) {
01191       log << MSG::DEBUG << "iTrack " << iTrack << endreq;
01192       aTrack = (*aRecMucTrackCol)[iTrack]; 
01193       //cout << "MucPos " << aTrack->GetMucPos() << " MucMomentum " << aTrack->GetMucMomentum() << endl;
01194 
01195       bool foundMcParticle = false;
01196       iter_mc = mcParticleCol->begin();
01197       for (;iter_mc != mcParticleCol->end(); iter_mc++) {
01198         if ((*iter_mc)->getTrackIndex() == aTrack->GetExtTrackID()) {
01199           log << MSG::DEBUG << " found McParticle index " << (*iter_mc)->getTrackIndex()
01200             << " equals to " << " Reconed track " << iTrack 
01201             << " index " << aTrack->GetExtTrackID() 
01202             << endreq;
01203           foundMcParticle = true;
01204           break;
01205         }
01206       }
01207       if (foundMcParticle == false) {
01208         nTracksMisFound++;
01209         log << MSG::DEBUG << " This Reconed track has no corresponding McParticle, where is it from? " << endreq; 
01210       }
01211 
01212       int nHitsFoundInMcParticle    = 0; 
01213       int nHitsMisFoundInMcParticle = 0; 
01214       vector< MucRecHit* > vecMucRecHit = aTrack->GetHits();
01215       log << MSG::DEBUG << " Reconed Track " << iTrack 
01216         << " index " << aTrack->GetExtTrackID() 
01217         << " contains " << aTrack->GetTotalHits() << " MucRecMcHit " << endreq;
01218       vector< MucRecHit* >::iterator iter_MucRecHit = vecMucRecHit.begin();
01219       for (; iter_MucRecHit != vecMucRecHit.end(); iter_MucRecHit++) {
01220         int part  = (*iter_MucRecHit)->Part();
01221         int seg   = (*iter_MucRecHit)->Seg();
01222         int gap   = (*iter_MucRecHit)->Gap();
01223         int strip = (*iter_MucRecHit)->Strip();
01224 
01225         log << MSG::DEBUG 
01226           << " Reconed track index " 
01227           << aTrack->GetExtTrackID() 
01228           << " contains " << " MucRecHit " 
01229           << " part "  << part
01230           << " seg "   << seg
01231           << " gap "   << gap
01232           << " strip " << strip;
01233 
01234         bool foundMcHit = false;
01235         vector< McPartToMucHitRel* > vecMucMcHit = assocMcMuc.getRelByFirst(*iter_mc);
01236         vector< McPartToMucHitRel* >::iterator iter_MucMcHit = vecMucMcHit.begin();
01237         for (; iter_MucMcHit != vecMucMcHit.end(); iter_MucMcHit++) {
01238           mucID = (*iter_MucMcHit)->getSecond()->identify();
01239           if (part  == MucID::barrel_ec(mucID) &&
01240               seg   == MucID::segment(mucID) &&
01241               gap   == MucID::layer(mucID) &&
01242               strip == MucID::channel(mucID)) {
01243             foundMcHit = true;
01244             break;
01245           }
01246         }
01247 
01248         if (foundMcHit == true) {
01249           nHitsFoundInMcParticle++;
01250           log << MSG::DEBUG << " This MucRecHit belongs to this Mc Particle " << endreq;
01251         }
01252         else { 
01253           nHitsMisFoundInMcParticle++;
01254           log << MSG::DEBUG << " This MucRecHit does not belong to this Reco track, it should not be attached " << endreq;
01255         }
01256       }
01257 
01258       //nHitsFound += nHitsFoundInMcParticle;
01259       nHitsMisFound  += nHitsMisFoundInMcParticle;
01260 
01261       log << MSG::DEBUG << " In " << vecMucRecHit.size() << " MucRecHit : "
01262         << nHitsFoundInMcParticle << " found in Mc Particle, " 
01263         << nHitsMisFoundInMcParticle << " not found in Mc Particle, mis attached " 
01264         << endreq;
01265     }
01266   }
01267 */
01268 /*
01269   log << MSG::INFO << " This event contains " << nTracksTotal << " Mc Particle, " 
01270     << nTracksFound << " tracks found, " 
01271     << nTracksLost  << " tracks lost, " 
01272     << nTracksMisFound << " tracks mis found " 
01273     << endreq;
01274   log << MSG::INFO << " This event contains " << nDigisTotal << " Muc Digis " << endreq;
01275   log << MSG::INFO << " This event contains " << nHitsTotal  << " MucMcHits, " 
01276     << nHitsFound << " mc hits found, "
01277     << nHitsLost  << " mc hits lost, " 
01278     << nHitsMisFound << " mc hits mis found "
01279     << endreq;
01280 */
01281   m_NDigisTotal        += nDigisTotal;
01282   m_NHitsTotal         += nHitsTotal;
01283   m_NHitsFoundTotal    += nHitsFound;
01284   m_NHitsLostTotal     += nHitsLost;
01285   m_NHitsMisFoundTotal += nHitsMisFound;
01286   //m_NHitsLost[nHitsLost]++;
01287 
01288   m_NTracksTotal          += nTracksTotal;
01289   m_NTracksRecoTotal      += nTracksFound;
01290   m_NTracksLostByMdcTotal += nTracksLost;
01291   m_NTracksLostByExtTotal += nTracksLostByExt;
01292   m_NTracksMdcGhostTotal  += nTracksMisFound;
01293   if (aRecMucTrackCol->size() > 0) 
01294   {
01295     RecMucTrack *aRecMucTrack = (*aRecMucTrackCol)[0]; 
01296     //m_posx = 0.0; m_posy = 0.0; m_posz = 0.0; 
01297     if(m_NtOutput>=1){
01298       m_posx = aRecMucTrack->getMucPos().x();
01299       m_posy = aRecMucTrack->getMucPos().y();
01300       m_posz = aRecMucTrack->getMucPos().z();
01301 
01302       m_posx_sigma = aRecMucTrack->getMucPosSigma().x();
01303       m_posy_sigma = aRecMucTrack->getMucPosSigma().y();
01304       m_posz_sigma = aRecMucTrack->getMucPosSigma().z();
01305     }
01306     //cout << m_posx << " " << m_posy << " " << m_posz << endl;
01307     //    m_tuple->write();
01308   }
01309   
01310   if(m_NtOutput>=1) m_tuple->write(); 
01311   return StatusCode::SUCCESS;
01312 }

StatusCode MucRecTrkExt::finalize (   ) 
 

StatusCode MucRecTrkExt::finalize (   ) 
  

01315                                   {
01316 
01317   MsgStream log(msgSvc(), name());
01318   log << MSG::INFO << "in finalize()" << endreq;
01319 
01320   log << MSG::INFO << " In " << m_NDigisTotal << " total MucDigis " << endreq;
01321   log << MSG::INFO << " In " << m_NHitsTotal  << " total MucMcHits " << endreq;
01322   log << MSG::INFO << m_NHitsFoundTotal    << " hits found in Reco tracks, " << endreq;
01323   log << MSG::INFO << m_NHitsLostTotal     << " hits lost (not found), of which " 
01324       << m_NHitsLostByMdcTotal << " lost because Mdc track not constructed " 
01325       << m_NHitsLostByExtTotal << " lost because Ext track not intersect with muc " << endreq;
01326   log << MSG::INFO << m_NHitsMisFoundTotal << " hits mis found (not belong to this track, but mis attached)" << endreq; 
01327   
01328   log << MSG::INFO << " In " << m_NTracksTotal << " total Mc tracks " << endreq;
01329   log << MSG::INFO << m_NTracksRecoTotal << " tracks found " << endreq;
01330   log << MSG::INFO << m_NTracksLostByMdcTotal << " tracks lost because Mdc track not constructed " << endreq;
01331   log << MSG::INFO << m_NTracksLostByExtTotal << " tracks lost because Ext track not intersect with muc " << endreq;
01332   log << MSG::INFO << m_NTracksMdcGhostTotal  << " tracks are Mdc ghost tracks " << endreq;
01333  
01334   /*
01335   for(int i = 0; i < 20; i++) log << MSG::DEBUG << "lost " << i << " hits track " << m_NHitsLost[i] << endreq;
01336   for(int i = 0; i < 9; i++)  log << MSG::DEBUG << "lost on gap " << i << " is " << m_NHitsLostInGap[i] << endreq;
01337   cout << m_TrackLostHit.size() << " track has hit lost, their event id : " << endl;
01338   for (int i = 0; i < m_TrackLostHit.size(); i++) {
01339     cout << m_TrackLostHit[i] << " : " << m_TrackLostHitNumber[i] << endl;
01340   }
01341   cout << m_TrackLostByMdc.size() << " tracks lost by Mdc, their event id : " << endl;
01342   for (int i = 0; i < m_TrackLostByMdc.size(); i++) {
01343     cout << m_TrackLostByMdc[i] << endl;
01344   }
01345   cout << m_TrackLostByExt.size() << " tracks lost by Ext no intersection with muc, their event id : " << endl;
01346   for (int i = 0; i < m_TrackLostByExt.size(); i++) {
01347     cout << m_TrackLostByExt[i] << endl;
01348   }
01349   */
01350 
01351   return StatusCode::SUCCESS;
01352 }

float MucRecTrkExt::getWindowSize (  Hep3Vector  mom,
int  part,
int  seg,
int  gap
) 
 

float MucRecTrkExt::getWindowSize (  Hep3Vector  mom,
int  part,
int  seg,
int  gap
) 
  

01482 {
01483   int mom_id;
01484   if(mom.mag()<0.6) mom_id = 0;
01485   else if(mom.mag()<0.8) mom_id = 1;
01486   else if(mom.mag()<1) mom_id = 2;
01487   else mom_id = 3;
01488 
01489   return kWindowSize_Mom_Gap[mom_id][gap]; 
01490 }

StatusCode MucRecTrkExt::initialize (   ) 
 

StatusCode MucRecTrkExt::initialize (   ) 
  

00066 {
00067   MsgStream log(msgSvc(), name());
00068   log << MSG::INFO << "in initialize()" << endreq;
00069 
00070   // Get the Particle Properties Service
00071   IPartPropSvc* p_PartPropSvc;
00072   static const bool CREATEIFNOTTHERE(true);
00073   StatusCode PartPropStatus = service("PartPropSvc", p_PartPropSvc, CREATEIFNOTTHERE);
00074   if (!PartPropStatus.isSuccess() || 0 == p_PartPropSvc) {
00075     log << MSG::WARNING << " Could not initialize Particle Properties Service" << endreq;
00076     return PartPropStatus;
00077   }
00078   m_particleTable = p_PartPropSvc->PDT();
00079   
00080   m_totalEvent          = 0;
00081   
00082   m_NDigisTotal         = 0;
00083   m_NHitsTotal          = 0;
00084   m_NHitsFoundTotal     = 0;
00085   m_NHitsLostTotal      = 0;
00086   m_NHitsMisFoundTotal  = 0;
00087   m_NHitsLostByMdcTotal = 0;  
00088   m_NHitsLostByExtTotal = 0;  
00089 
00090   m_NTracksTotal     = 0;
00091   m_NTracksRecoTotal = 0;
00092   m_NTracksLostByMdcTotal = 0;
00093   m_NTracksLostByExtTotal = 0;
00094   m_NTracksMdcGhostTotal  = 0;
00095 
00096   for(int i = 0; i < 20; i++) m_NHitsLost.push_back(0);
00097   for(int i = 0; i < 10; i++) m_NHitsLostInGap.push_back(0);
00098 
00099   IMucGeomSvc* mucGeomSvc;
00100   StatusCode sc = service("MucGeomSvc", mucGeomSvc);
00101   if (sc ==  StatusCode::SUCCESS) {              
00102     //cout <<"dump"<<endl;
00103    mucGeomSvc->Dump();
00104    //cout<<"Hi, event routine is running"<<endl;
00105   } else {
00106     return StatusCode::FAILURE;
00107   }
00108   m_MucRecHitContainer = 0;
00109 
00110   //--------------book ntuple------------------
00111   //  NTuplePtr nt1(ntupleSvc(), "MyTuples/1");
00112   if(m_NtOutput>=1){
00113     NTuplePtr nt1(ntupleSvc(), "FILE401/T");
00114 
00115     if ( nt1 ) { m_tuple = nt1;}
00116     else {
00117       //    m_tuple = ntupleSvc()->book ("MyTuples/1", CLID_RowWiseTuple, "MdcTrkRecon N-Tuple");
00118       m_tuple = ntupleSvc()->book ("FILE401/T", CLID_RowWiseTuple, "MucTrkRecon N-Tuple");
00119       if ( m_tuple )    {
00120         sc = m_tuple->addItem ("posx",    m_posx);
00121         sc = m_tuple->addItem ("posy",    m_posy);
00122         sc = m_tuple->addItem ("posz",    m_posz);
00123         sc = m_tuple->addItem ("posx_ext",    m_posx_ext);
00124         sc = m_tuple->addItem ("posy_ext",    m_posy_ext);
00125         sc = m_tuple->addItem ("posz_ext",    m_posz_ext);
00126         sc = m_tuple->addItem ("posxsigma",   m_posx_sigma);
00127         sc = m_tuple->addItem ("posysigma",   m_posy_sigma);
00128         sc = m_tuple->addItem ("poszsigma",   m_posz_sigma);
00129         sc = m_tuple->addItem ("Depth",   m_depth);
00130         sc = m_tuple->addItem ("Distance",m_Distance);
00131         sc = m_tuple->addItem ("MaxHits", m_MaxHits);
00132         sc = m_tuple->addItem ("Chi2",    m_Chi2);
00133         sc = m_tuple->addItem ("Dist_x",  m_Dist_x);
00134         sc = m_tuple->addItem ("Dist_y",  m_Dist_y);
00135         sc = m_tuple->addItem ("Dist_z",  m_Dist_z);
00136         sc = m_tuple->addItem ("px_mc",   m_px_mc);
00137         sc = m_tuple->addItem ("py_mc",   m_py_mc);
00138         sc = m_tuple->addItem ("pz_mc",   m_pz_mc);
00139         sc = m_tuple->addItem ("emctrack",m_emctrack);
00140         sc = m_tuple->addItem ("muchasdigi",m_muc_digi);
00141 
00142         sc = m_tuple->addItem ("part",  m_part);
00143         sc = m_tuple->addItem ("seg",   m_seg);
00144         sc = m_tuple->addItem ("gap",   m_gap);
00145         sc = m_tuple->addItem ("strip", m_strip);
00146         sc = m_tuple->addItem ("diff",  m_diff);
00147         sc = m_tuple->addItem ("distance",m_distance);
00148         sc = m_tuple->addItem ("run",   m_run);
00149         sc = m_tuple->addItem ("event", m_event);
00150 
00151       }
00152       else    {   // did not manage to book the N tuple....
00153         log << MSG::ERROR << "    Cannot book N-tuple:" << long(m_tuple) << endmsg;
00154         //return StatusCode::FAILURE;
00155       }
00156     }
00157   }
00158   //--------------end of book ntuple------------------
00159 
00160   return StatusCode::SUCCESS;
00161 }

void MucRecTrkExt::TrackFinding (  RecMucTrack *  aTrack  ) 
 

void MucRecTrkExt::TrackFinding (  RecMucTrack *  aTrack  ) 
  

01356 {
01357   MsgStream log(msgSvc(), name());
01358 
01359   Hep3Vector currentPos = aTrack->GetCurrentPos();
01360   Hep3Vector currentDir = aTrack->GetCurrentDir();
01361 //   if(currentPos.mag() < kMinor) {
01362 //     log << MSG::WARNING << "No MUC intersection in track " << endreq;
01363 //     continue;
01364 //   }
01365    
01366   int firstHitFound[2] = { 0,  0}; // Has the fist position in this orient determined? if so, could CorrectDirection()
01367   int firstHitGap[2]   = {-1, -1}; // When the fist position in this orient determined, the gap it is on
01368   for(int partSeq = 0; partSeq < (int)MucID::getPartNum(); partSeq++) 
01369   {
01370     int iPart = kPartSeq[partSeq];
01371     for(int iGap = 0; iGap < (int)MucID::getGapNum(iPart); iGap++) 
01372     {
01373       int seg = -1;
01374       int orient = MucGeoGeneral::Instance()->GetGap(iPart, 0, iGap)->Orient();;
01375       float xInsct, yInsct, zInsct;
01376       aTrack->Project(iPart, iGap, xInsct, yInsct, zInsct, seg);
01377       log << MSG::INFO << "part " << iPart << " gap " << iGap << " x " << xInsct << " y " << yInsct << " z " << zInsct << " seg " << seg << endreq; 
01378       //cout<<"project: gap="<<iGap<<"  seg="<<seg<<"  "<<xInsct<<" "<<yInsct<<" "<<zInsct<<endl;
01379 
01380       if(seg == -1) {
01381         //log << MSG::DEBUG << "no intersection with part " << iPart
01382         //  << " gap " << iGap << " in this track !" << endl;
01383         continue;
01384       }
01385       aTrack->SetCurrentInsct(xInsct, yInsct, zInsct);
01386       
01387       for(int iDeltaSeg = 0; iDeltaSeg < kNSegSearch; iDeltaSeg++) 
01388       {
01389         int iSeg = seg + kDeltaSeg[iDeltaSeg];  // also find on neighbor seg;
01390         if(iSeg < 0) iSeg += MucID::getSegNum(iPart);
01391         if(iSeg > (int)MucID::getSegNum(iPart) - 1) iSeg -= MucID::getSegNum(iPart);
01392 
01393         //float window = kWindowSize[iPart][iGap];
01394         Hep3Vector mom_mdc = aTrack->getMdcMomentum();
01395         float window = getWindowSize(mom_mdc, iPart, iSeg, iGap);
01396 
01397         if (firstHitFound[orient] != 1) window *= kFirstHitWindowRatio;   // if no hits have been found on this orient, expand the window.
01398 
01399         for (int iHit = 0; iHit < m_MucRecHitContainer->GetGapHitCount(iPart, iSeg, iGap); iHit++) 
01400         {
01401           log << MSG::DEBUG << "iSeg " << iSeg << " iHit " << iHit << endreq;
01402           MucRecHit* pHit = m_MucRecHitContainer->GetHit(iPart, iSeg, iGap, iHit);
01403           //cout<< "strip " << pHit->Strip() << " center " << pHit->GetCenterPos() << endl;
01404 
01405           if (!pHit) {
01406             log << MSG::WARNING << "MucRecTrkExt: null pointer to pHit" << endreq;
01407             continue;
01408           }
01409           else 
01410           {
01411             // Get the displacement of hit pHit to intersection
01412             float dX  = aTrack->GetHitDistance2(pHit);
01413             log << MSG::DEBUG << "distance = " << setw(8) << dX << " size " << setw(4) << window << endreq;
01414 
01415             //cout<<"dX= "<<dX<<"  window="<<window<<endl;
01416             if (fabs(dX) < window) 
01417             {
01418               // Attach the hit if it exists
01419               //cout << "meet window " << pHit->GetSoftID() << endl;
01420               //****************if this if emc track, we abondon used hit in mdc*******************
01421               //if(m_emcrec == 1 )
01422               if(aTrack->GetRecMode() == 0) {
01423                 pHit->SetHitMode(1); //mdc ext
01424                 aTrack->AttachHit(pHit);
01425                 //cout<<"in MucRecTrkExt: trackmode==0 so mdc ext "<<iPart<<" "<<iSeg<<" "<<iGap<<" "<<iHit<<endl;
01426               }
01427               if(aTrack->GetRecMode() == 1) {
01428                 //cout<<"in MucRecTrkExt: HitMode = "<<pHit->GetHitMode()<<" part:  "<<iPart<<" "<<iSeg<<" "<<iGap<<" "<<iHit<<endl;
01429                 if(pHit->GetHitMode()!=1) {
01430                   aTrack->AttachHit(pHit);  //this hit has not been used by mdc ext
01431                   pHit->SetHitMode(2); //emc ext
01432                 }
01433               }
01434               if(aTrack->GetRecMode() == 2) {
01435                 //cout<<"in MucRecTrkExt: HitMode = "<<pHit->GetHitMode()<<" part:  "<<iPart<<" "<<iSeg<<" "<<iGap<<" "<<iHit<<endl;
01436                 if(pHit->GetHitMode()==-1) {
01437                   aTrack->AttachHit(pHit);  //this hit has not been used by mdc ext
01438                   pHit->SetHitMode(2); //emc ext
01439                 }
01440               }
01441 
01442               if (firstHitGap[orient] == -1) firstHitGap[orient] = iGap;
01443               firstHitFound[orient] = 1;
01444               //cout << "You could correct directon in orient " << orient << endl;
01445               //cout<< " part " << iPart << " seg " << iSeg << " gap " << iGap 
01446               //          << " strip " << setw(2) << pHit->Strip() << " attatched" << endl;
01447               log << MSG::DEBUG << " part " << iPart << " seg " << iSeg << " gap " << iGap 
01448                 << " strip " << setw(2) << pHit->Strip() << " attatched" << endreq;
01449               log << MSG::DEBUG << "current total hits " << aTrack->GetTotalHits() << endreq;
01450             }
01451             else 
01452             {
01453               m_NHitsLostInGap[iGap]++;
01454               //log << MSG::DEBUG cout << " part " << iPart << " seg " << iSeg << " gap " << iGap 
01455               //       << " strip " << setw(2) << pHit->GetSoftID().GetStrip() 
01456               //     << " not attached !" << endreq;
01457             }
01458           } // end pHit else
01459         } // end iHit
01460         
01461         aTrack->CalculateInsct(iPart, iSeg, iGap);
01462       } // end DeltaSeg
01463       
01464       // When correct dir in the orient is permitted and this gap is not the gap first hit locates. 
01465       if (firstHitFound[orient] && firstHitGap[orient] != iGap) aTrack->CorrectDir();
01466       aTrack->CorrectPos();
01467       //cout << "Current Intersection " << aTrack->GetCurrentInsct() << endl;
01468       //cout << "Current Direction    " << aTrack->GetCurrentDir() << endl;
01469       aTrack->AttachInsct(aTrack->GetCurrentInsct());
01470     } // end iGap
01471   } // end iPart
01472   
01473   aTrack->LineFit(m_fittingMethod);
01474   aTrack->ComputeTrackInfo(m_fittingMethod);
01475 
01476   //cout<<" in TrackFinding: depth= "<<aTrack->GetDepth()<<endl;
01477 } // End TrackFinding()

Member Data Documentation

bool MucRecTrkExt::m_Blind [private]
 

NTuple::Item<double> MucRecTrkExt::m_Chi2 [private]
 

NTuple::Item<double> MucRecTrkExt::m_Chi2 [private]
 

int MucRecTrkExt::m_CompareWithMcHit [private]
 

std::string MucRecTrkExt::m_configFile [private]
 

NTuple::Item<double> MucRecTrkExt::m_depth [private]
 

NTuple::Item<double> MucRecTrkExt::m_depth [private]
 

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

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

NTuple::Item<double> MucRecTrkExt::m_Dist_x [private]
 

NTuple::Item<double> MucRecTrkExt::m_Dist_x [private]
 

NTuple::Item<double> MucRecTrkExt::m_Dist_y [private]
 

NTuple::Item<double> MucRecTrkExt::m_Dist_y [private]
 

NTuple::Item<double> MucRecTrkExt::m_Dist_z [private]
 

NTuple::Item<double> MucRecTrkExt::m_Dist_z [private]
 

NTuple::Item<double> MucRecTrkExt::m_distance [private]
 

NTuple::Item<double> MucRecTrkExt::m_Distance [private]
 

NTuple::Item<double> MucRecTrkExt::m_distance [private]
 

NTuple::Item<double> MucRecTrkExt::m_Distance [private]
 

int MucRecTrkExt::m_emcrec [private]
 

int MucRecTrkExt::m_EmcShowerSeedMode [private]
 

NTuple::Item<double> MucRecTrkExt::m_emctrack [private]
 

NTuple::Item<double> MucRecTrkExt::m_emctrack [private]
 

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

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

int MucRecTrkExt::m_ExtTrackSeedMode [private]
 

int MucRecTrkExt::m_fittingMethod [private]
 

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

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

NTuple::Item<double> MucRecTrkExt::m_MaxHits [private]
 

NTuple::Item<double> MucRecTrkExt::m_MaxHits [private]
 

bool MucRecTrkExt::m_MsOutput [private]
 

NTuple::Item<double> MucRecTrkExt::m_muc_digi [private]
 

NTuple::Item<double> MucRecTrkExt::m_muc_digi [private]
 

int MucRecTrkExt::m_MucHitSeedMode [private]
 

MucRecHitContainer* MucRecTrkExt::m_MucRecHitContainer [private]
 

MucRecHitContainer* MucRecTrkExt::m_MucRecHitContainer [private]
 

int MucRecTrkExt::m_NDigisTotal [private]
 

int MucRecTrkExt::m_NHitsFoundTotal [private]
 

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

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

int MucRecTrkExt::m_NHitsLostByExtTotal [private]
 

int MucRecTrkExt::m_NHitsLostByMdcTotal [private]
 

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

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

int MucRecTrkExt::m_NHitsLostTotal [private]
 

int MucRecTrkExt::m_NHitsMisFoundTotal [private]
 

int MucRecTrkExt::m_NHitsTotal [private]
 

int MucRecTrkExt::m_NtOutput [private]
 

int MucRecTrkExt::m_NTracksLostByExtTotal [private]
 

int MucRecTrkExt::m_NTracksLostByMdcTotal [private]
 

int MucRecTrkExt::m_NTracksMdcGhostTotal [private]
 

int MucRecTrkExt::m_NTracksRecoTotal [private]
 

int MucRecTrkExt::m_NTracksTotal [private]
 

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

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

HepPDT::ParticleDataTable* MucRecTrkExt::m_particleTable [private]
 

HepPDT::ParticleDataTable* MucRecTrkExt::m_particleTable [private]
 

NTuple::Item<double> MucRecTrkExt::m_posx [private]
 

NTuple::Item<double> MucRecTrkExt::m_posx [private]
 

NTuple::Item<double> MucRecTrkExt::m_posx_ext [private]
 

NTuple::Item<double> MucRecTrkExt::m_posx_ext [private]
 

NTuple::Item<double> MucRecTrkExt::m_posx_sigma [private]
 

NTuple::Item<double> MucRecTrkExt::m_posx_sigma [private]
 

NTuple::Item<double> MucRecTrkExt::m_posy [private]
 

NTuple::Item<double> MucRecTrkExt::m_posy [private]
 

NTuple::Item<double> MucRecTrkExt::m_posy_ext [private]
 

NTuple::Item<double> MucRecTrkExt::m_posy_ext [private]
 

NTuple::Item<double> MucRecTrkExt::m_posy_sigma [private]
 

NTuple::Item<double> MucRecTrkExt::m_posy_sigma [private]
 

NTuple::Item<double> MucRecTrkExt::m_posz [private]
 

NTuple::Item<double> MucRecTrkExt::m_posz [private]
 

NTuple::Item<double> MucRecTrkExt::m_posz_ext [private]
 

NTuple::Item<double> MucRecTrkExt::m_posz_ext [private]
 

NTuple::Item<double> MucRecTrkExt::m_posz_sigma [private]
 

NTuple::Item<double> MucRecTrkExt::m_posz_sigma [private]
 

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

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

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

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

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

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

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

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

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

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

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

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

long MucRecTrkExt::m_totalEvent [private]
 

std::vector<int> MucRecTrkExt::m_TrackLostByExt [private]
 

std::vector<int> MucRecTrkExt::m_TrackLostByExt [private]
 

std::vector<int> MucRecTrkExt::m_TrackLostByMdc [private]
 

std::vector<int> MucRecTrkExt::m_TrackLostByMdc [private]
 

std::vector<int> MucRecTrkExt::m_TrackLostHit [private]
 

std::vector<int> MucRecTrkExt::m_TrackLostHit [private]
 

std::vector<int> MucRecTrkExt::m_TrackLostHitNumber [private]
 

std::vector<int> MucRecTrkExt::m_TrackLostHitNumber [private]
 

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

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

The documentation for this class was generated from the following files:
Generated on Wed Feb 2 16:36:25 2011 for BOSS6.5.5 by  doxygen 1.3.9.1