BOSS6.5.5: MdcxCosmicSewer Class Reference

00056                                                                                 :
00057   Algorithm(name, pSvcLocator) 
00058 {
00059 
00060   declareProperty("debug",          m_debug = false);
00061   declareProperty("hist",           m_hist  = false);
00062 
00063   declareProperty("doSag",          m_doSag = false);
00064   
00065   declareProperty("cosmicSewPar",   m_cosmicSewPar);
00066   declareProperty("cosmicSewSkip",  m_cosmicSewSkip=false);
00067 
00068   declareProperty("countPropTime",  m_countPropTime = true);
00069 }

00074                                   {
00075   delete m_bfield;
00076 }

00078                                     {  
00079   //Initailize MdcDetector
00080   m_gm = MdcDetector::instance(m_doSag); 
00081   if(NULL == m_gm) return StatusCode::FAILURE;
00082 
00083   //Initailize MdcCheckUtil
00084   m_mdcCheckUtil = new MdcCheckUtil();
00085 
00086   return StatusCode::SUCCESS;
00087 }

00471                                          {
00472   MsgStream log(msgSvc(), name());
00473   StatusCode sc;
00474   //-------------------------------------------------------
00475   //clear TDS
00476   //-------------------------------------------------------
00477   IDataManagerSvc *dataManSvc;
00478   dataManSvc = dynamic_cast<IDataManagerSvc*> (eventSvc()); 
00479   DataObject *aTrackCol;
00480   eventSvc()->findObject("/Event/Recon/RecMdcTrackCol",aTrackCol);
00481   if(aTrackCol != NULL) {
00482     dataManSvc->clearSubTree("/Event/Recon/RecMdcTrackCol");
00483     eventSvc()->unregisterObject("/Event/Recon/RecMdcTrackCol");
00484   }
00485   DataObject *aRecHitCol;
00486   eventSvc()->findObject("/Event/Recon/RecMdcHitCol",aRecHitCol);
00487   if(aRecHitCol != NULL) {
00488     dataManSvc->clearSubTree("/Event/Recon/RecMdcHitCol");
00489     eventSvc()->unregisterObject("/Event/Recon/RecMdcHitCol");
00490   }
00491   SmartDataPtr<RecMdcTrackCol> trackListTemp(eventSvc(),EventModel::Recon::RecMdcTrackCol);
00492   if (!trackListTemp) {
00493     RecMdcTrackCol *trackListTemp= new RecMdcTrackCol;
00494     sc =  eventSvc()->registerObject(EventModel::Recon::RecMdcTrackCol, trackListTemp);
00495     if(!sc.isSuccess()) {
00496       log << MSG::FATAL << " Could not register RecMdcTrack collection" <<endreq;
00497     }
00498   }
00499   SmartDataPtr<RecMdcHitCol> hitListTemp(eventSvc(),EventModel::Recon::RecMdcHitCol);
00500   if (!hitListTemp) {
00501     RecMdcHitCol *hitListTemp= new RecMdcHitCol;
00502     sc =  eventSvc()->registerObject(EventModel::Recon::RecMdcHitCol, hitListTemp);
00503     if(!sc.isSuccess()) {
00504       log << MSG::FATAL << " Could not register RecMdcHit collection" <<endreq;
00505     }
00506   }
00507 }

00144                                     {
00145   MsgStream log(msgSvc(), name());
00146   log << MSG::INFO << "in execute()" << endreq; 
00147 
00148 
00149   // Get event No.
00150   SmartDataPtr<Event::EventHeader> evtHead(eventSvc(),"/Event/EventHeader");
00151   if (!evtHead) {
00152     log << MSG::FATAL<< "Could not retrieve event header" << endreq;
00153     return StatusCode::FAILURE;
00154   } 
00155   long t_evtNo = evtHead->eventNumber();
00156   if (m_debug) std::cout << "sew evt: " << t_evtNo << std::endl;
00157 
00158 
00159   // Get bunch time t0 (ns)
00160   m_bunchT0 = -999.;
00161   SmartDataPtr<RecEsTimeCol> aevtimeCol(eventSvc(),"/Event/Recon/RecEsTimeCol");
00162   if (!aevtimeCol || aevtimeCol->size()==0) {
00163     log << MSG::WARNING<< "Could not find RecEsTimeCol"<< endreq;
00164     return StatusCode::SUCCESS;
00165   }
00166   RecEsTimeCol::iterator iter_evt = aevtimeCol->begin();
00167   for(; iter_evt!=aevtimeCol->end(); iter_evt++){
00168     m_bunchT0 =  (*iter_evt)->getTest();
00169   }
00170 
00171   //read event data 
00172   SmartDataPtr<RecMdcTrackCol> recMdcTrackCol(eventSvc(),EventModel::Recon::RecMdcTrackCol);
00173   SmartDataPtr<RecMdcHitCol> recMdcHitCol(eventSvc(),EventModel::Recon::RecMdcHitCol);
00174   if (!recMdcTrackCol || ! recMdcHitCol) return StatusCode::SUCCESS; 
00175   if (2!=recMdcTrackCol->size()){
00176     clearRecMdcTrackHit();
00177     return StatusCode::SUCCESS;
00178   }
00179   if(m_debug)std::cout<<__FILE__<<" nTk==2 begin sewing"<<  std::endl;
00180 
00181 
00182   //get best track
00183   double dParam[5]={0.,0.,0.,0.,0.};
00184   float bprob =0.;
00185   float chisum =0.; 
00186   RecMdcTrackCol::iterator besthel;
00187 
00188   RecMdcTrackCol::iterator it = recMdcTrackCol->begin();
00189   for (int iTk=0; it != recMdcTrackCol->end(); it++,iTk++) {
00190     RecMdcTrack* tk = *it;
00191     double Bz = m_bfield->bFieldNominal();
00192     double d0 = -tk->helix(0);    //cm
00193     double phi0 = tk->helix(1)+ CLHEP::halfpi;
00194     double omega = Bz*tk->helix(2)/-333.567;
00195     double z0 = tk->helix(3);    //cm
00196     double tanl = tk->helix(4);    //cm
00197 
00198     if(phi0>Constants::pi)phi0-=Constants::twoPi;
00199     if(phi0<-Constants::pi)phi0+=Constants::twoPi;
00200 
00201     if(m_debug)std::cout<<__FILE__<<" "<<__LINE__<< "tk"<<iTk
00202       <<"("<<d0<<","<<phi0<<"," <<","<<omega<<","<<z0<<","<<tanl<<")"<<std::endl;
00203 
00204     if(iTk==0){//store param1
00205       dParam[0]=d0;
00206       dParam[1]=phi0;
00207       dParam[2]=omega;
00208       dParam[3]=z0;
00209       dParam[4]=tanl;
00210     }else{//calc. diff between param1 and param2
00211       dParam[0] += d0;
00212       dParam[1] -= (phi0+Constants::pi);
00213       dParam[2] -= omega;
00214       dParam[3] -= z0;
00215       dParam[4] += tanl;
00216     }
00217 
00218     if(phi0<0) besthel = it;
00219     float bchisq = tk->chi2();
00220     int   bndof  = tk->ndof();
00221     bprob=Mdcxprobab(bndof,bchisq);
00222     chisum += bchisq;
00223     //if (bprob > bestprob){
00224     //bestprob = bprob; 
00225     //besthel = it;
00226     //}
00227   }//end for recMdcTrackCol
00228 
00229   TrkHelixMaker helixfactory;
00230 
00231   if(m_debug){
00232     std::cout<<__FILE__<<" param diff: " << "\n D0 " << dParam[0]
00233       << "\n Phi0 " << dParam[1] << "\n Omega " << dParam[2]
00234       << "\n Z0 " << dParam[3] << "\n Tanl " << dParam[4] << std::endl;
00235   }
00236 
00237   if(m_hist){
00238     m_csmcD0= dParam[0];
00239     m_csmcPhi0=dParam[1];
00240     m_csmcOmega=dParam[2];
00241     m_csmcZ0=dParam[3];
00242     m_csmcTanl=dParam[4];
00243     m_xtupleCsmcSew->write();
00244   }
00245 
00246   //get track 5 parameters
00247   if( (fabs(dParam[0])>m_cosmicSewPar[0]) ||
00248       (fabs(dParam[1])>m_cosmicSewPar[1]) || 
00249       (fabs(dParam[2])>m_cosmicSewPar[2]) ||
00250       (fabs(dParam[3])>m_cosmicSewPar[3]) ||
00251       (fabs(dParam[4])>m_cosmicSewPar[4]) ){
00252     if(m_debug)std::cout<<__FILE__<<" 2 trk param not satisfied skip!"<<std::endl;
00253     clearRecMdcTrackHit();
00254     return StatusCode::SUCCESS;
00255   }
00256 
00257   //got 2 sew-able trks, get helix param
00258   RecMdcTrack* tk = *besthel;
00259   double Bz = m_bfield->bFieldNominal();
00260   double d0 = -tk->helix(0);    //cm
00261   double phi0 = tk->helix(1)+ CLHEP::halfpi;
00262   double omega = Bz*tk->helix(2)/-333.567;
00263   double z0 = tk->helix(3);    //cm
00264   double tanl = tk->helix(4);    //cm
00265 
00266   TrkExchangePar tt(d0,phi0,omega,z0,tanl);
00267   if(m_debug)std::cout<<__FILE__<<" besthel:"<<d0<<" "<<phi0
00268     <<" "<<omega<<" "<<z0<<" "<<tanl<<  std::endl;
00269 
00270 
00271   //new track
00272   TrkRecoTrk* newTrack = helixfactory.makeTrack(tt, chisum, *m_context, m_bunchT0*1.e-9);
00273   helixfactory.setFlipAndDrop(*newTrack, true, true);
00274 
00275 
00276   // combine hit list
00277   TrkHitList* qhits = newTrack->hits();
00278   HitRefVec skippedHits; 
00279   it = recMdcTrackCol->begin();
00280   for (int iTk=0; it != recMdcTrackCol->end(); it++,iTk++) {
00281     RecMdcTrack* tk = *it;
00282     HitRefVec hits = tk->getVecHits();
00283     MdcxHitsToHots(qhits,hits,skippedHits);
00284   }
00285 
00286 
00287   //------------------------------------
00288   //do fit
00289   //------------------------------------
00290   TrkErrCode err=qhits->fit();
00291   const TrkFit* newFit = newTrack->fitResult(); 
00292 
00293 
00294   if (!newFit) {
00295     // Fit bad
00296     if(m_debug)std::cout<<__FILE__<<" sew fit failed!!!"<<  std::endl;
00297     delete newTrack;
00298 
00299   } else {
00300     //------------------------------------
00301     // Fit good
00302     //------------------------------------
00303     helixfactory.setFlipAndDrop(*newTrack, false, false);
00304     if(m_debug){
00305       std::cout<<__FILE__<<" sew cosmic fit good "<<  std::endl;
00306       newTrack->print(std::cout);
00307     }
00308 
00309     //calc. doca of skipped hits
00310     HitRefVec::iterator itHitSkipped = skippedHits.begin();   
00311     for (;itHitSkipped!=skippedHits.end();++itHitSkipped){
00312       RecMdcHit* h = *itHitSkipped;
00313       Identifier id(h->getMdcId());
00314       int layer = MdcID::layer(id);
00315       int wire  = MdcID::wire(id);
00316       double fltLen = 0.;
00317       HepVector helix = newFit->helix(fltLen).params();
00318       HepSymMatrix err = newFit->helix(fltLen).covariance();
00319       double doca = m_mdcCheckUtil->docaPatPar(layer,wire,helix,err);
00320       h->setDoca(doca);
00321     }
00322 
00323     //-------------------------------------------------------
00324     //Store TDS
00325     //-------------------------------------------------------
00326     IDataManagerSvc *dataManSvc;
00327     dataManSvc = dynamic_cast<IDataManagerSvc*> (eventSvc()); 
00328     DataObject *aTrackCol;
00329     eventSvc()->findObject("/Event/Recon/RecMdcTrackCol",aTrackCol);
00330     if(aTrackCol != NULL) {
00331       dataManSvc->clearSubTree("/Event/Recon/RecMdcTrackCol");
00332       eventSvc()->unregisterObject("/Event/Recon/RecMdcTrackCol");
00333     }
00334     DataObject *aRecHitCol;
00335     eventSvc()->findObject("/Event/Recon/RecMdcHitCol",aRecHitCol);
00336     if(aRecHitCol != NULL) {
00337       dataManSvc->clearSubTree("/Event/Recon/RecMdcHitCol");
00338       eventSvc()->unregisterObject("/Event/Recon/RecMdcHitCol");
00339     }
00340     store(newTrack, skippedHits);//newTrack have been deleted
00341   }// end if newFit
00342 
00343   return StatusCode::SUCCESS;
00344 }

00346                                      {
00347   MsgStream log(msgSvc(), name());
00348   log << MSG::INFO << "in finalize()" << endreq;
00349   return StatusCode::SUCCESS;
00350 }

00092                                       {  
00093   MsgStream log(msgSvc(), name());
00094   StatusCode sc;
00095   log << MSG::INFO << "in initialize()" << endreq;      
00096 
00097   //Initailize magnetic filed 
00098   sc = service ("MagneticFieldSvc",m_pIMF);
00099   if(sc!=StatusCode::SUCCESS) {
00100     log << MSG::ERROR << "Unable to open Magnetic field service"<<endreq;
00101   }
00102   m_bfield = new BField(m_pIMF);
00103   log << MSG::INFO << "field z = "<<m_bfield->bFieldNominal()<< endreq;
00104   m_context = new TrkContextEv(m_bfield);
00105 
00106 
00107   // Get MdcCalibFunSvc
00108   IMdcCalibFunSvc* imdcCalibSvc; 
00109   sc = service ("MdcCalibFunSvc", imdcCalibSvc);
00110   m_mdcCalibFunSvc = dynamic_cast<MdcCalibFunSvc*>(imdcCalibSvc);
00111   if ( sc.isFailure() ){
00112     log << MSG::FATAL << "Could not load MdcCalibFunSvc!" << endreq;
00113     return StatusCode::FAILURE;
00114   }
00115 
00116   // Get RawDataProviderSvc
00117   IRawDataProviderSvc* iRawDataProvider; 
00118   sc = service ("RawDataProviderSvc", iRawDataProvider);
00119   if ( sc.isFailure() ){
00120     log << MSG::FATAL << "Could not load RawDataProviderSvc!" << endreq;
00121     return StatusCode::FAILURE;
00122   }
00123   m_rawDataProviderSvc = dynamic_cast<RawDataProviderSvc*>(iRawDataProvider);
00124 
00125 
00126   //book ntuple
00127   NTuplePtr ntCsmcSew(ntupleSvc(), "FILE121/csmc");
00128   if ( ntCsmcSew ) { m_xtupleCsmcSew = ntCsmcSew;}
00129   else {
00130     m_xtupleCsmcSew = ntupleSvc()->book ("FILE121/csmc", CLID_ColumnWiseTuple, "MdcxReco reconsturction results");
00131     if ( m_xtupleCsmcSew ) {
00132       sc = m_xtupleCsmcSew->addItem ("dD0",      m_csmcD0);
00133       sc = m_xtupleCsmcSew->addItem ("dPhi0",    m_csmcPhi0);
00134       sc = m_xtupleCsmcSew->addItem ("dZ0",      m_csmcZ0);
00135       sc = m_xtupleCsmcSew->addItem ("dOmega",   m_csmcOmega);
00136       sc = m_xtupleCsmcSew->addItem ("dPt",      m_csmcPt);
00137       sc = m_xtupleCsmcSew->addItem ("dTanl",    m_csmcTanl);
00138     }
00139   }//end of book
00140 
00141   return StatusCode::SUCCESS;
00142 }

00352                                                                                                             {
00353 
00354   //get MdcDigi pointer
00355   uint32_t getDigiFlag = 0;
00356   MdcDigiVec mdcDigiVec = m_rawDataProviderSvc->getMdcDigiVec(getDigiFlag);
00357   if (0 == mdcDigiVec.size()){
00358     std::cout << " No hits in MdcDigiVec" << std::endl;
00359     return; 
00360   }
00361   const MdcDigi* m_digiMap[43][288];
00362   MdcDigiVec::iterator iter = mdcDigiVec.begin();
00363   for (; iter != mdcDigiVec.end(); iter++ ) {
00364     const MdcDigi* aDigi = *iter;
00365     const Identifier id= aDigi->identify();
00366     int layer = MdcID::layer(id);
00367     int wire = MdcID::wire(id);
00368     m_digiMap[layer][wire] = aDigi;
00369   }
00370 
00371 
00372   // new MdcRecoHitOnTracks
00373   HitRefVec::iterator itHit = recMdcHits.begin();   
00374   for (;itHit!=recMdcHits.end();++itHit){
00375     RecMdcHit* h = *itHit;
00376     Identifier id(h->getMdcId());
00377     int layer = MdcID::layer(id);
00378     int wire  = MdcID::wire(id);
00379     int flagLR= h->getFlagLR();
00380     int ambig = 2;
00381     if ( flagLR == 0 ){
00382       ambig = 1; //left driftDist <0
00383     }else if( flagLR == 1){
00384       ambig = -1; //right driftDist >0
00385     }else if( flagLR == 2){
00386       ambig = 0; //don't know 
00387     }
00388     MdcHit *thehit = new MdcHit(m_digiMap[layer][wire], m_gm); 
00389     thehit->setCalibSvc(m_mdcCalibFunSvc);
00390     thehit->setCountPropTime(m_countPropTime);
00391     MdcRecoHitOnTrack* newhot= new MdcRecoHitOnTrack(*thehit, ambig, m_bunchT0);//m_bunchT0 nano second here
00392     //MdcHitOnTrack* newhot = &temp;
00393     double fltLen = h->getFltLen();
00394     newhot->setFltLen(fltLen);
00395     if(m_cosmicSewSkip){ 
00396       RecMdcHit* newSkippedHit = new RecMdcHit();
00397       newSkippedHit->setId(h->getId());
00398       newSkippedHit->setTrkId(-1);
00399       newSkippedHit->setDriftDistLeft(h->getDriftDistLeft());
00400       newSkippedHit->setDriftDistRight(h->getDriftDistRight());
00401       newSkippedHit->setErrDriftDistLeft(h->getErrDriftDistLeft());
00402       newSkippedHit->setErrDriftDistRight(h->getErrDriftDistRight());
00403       newSkippedHit->setChisqAdd(h->getChisqAdd());
00404       newSkippedHit->setFlagLR(h->getFlagLR());
00405       newSkippedHit->setStat(h->getStat());
00406       newSkippedHit->setMdcId(h->getMdcId());
00407       newSkippedHit->setTdc(h->getTdc());
00408       newSkippedHit->setAdc(h->getAdc());
00409       newSkippedHit->setDriftT(h->getDriftT());
00410       newSkippedHit->setDoca(999.);
00411       newSkippedHit->setEntra(h->getEntra());
00412       newSkippedHit->setZhit(h->getZhit());
00413       newSkippedHit->setFltLen(h->getFltLen());
00414 
00415       if(layer<20) skippedHits.push_back(newSkippedHit);
00416       else m_trkHitList->appendHot(newhot);
00417     }else{
00418       m_trkHitList->appendHot(newhot);
00419     }
00420   }
00421 }

00423                                                                       {
00424   StatusCode sc;
00425   MsgStream log(msgSvc(), name());
00426 
00427   assert (aTrack != NULL);
00428 
00429   RecMdcTrackCol *trackList = new RecMdcTrackCol;
00430   sc =  eventSvc()->registerObject(EventModel::Recon::RecMdcTrackCol, trackList);
00431   if(!sc.isSuccess()) {
00432     log << MSG::FATAL << " Could not register RecMdcTrack collection" <<endreq;
00433   }
00434   RecMdcHitCol *hitList = new RecMdcHitCol;
00435   sc =  eventSvc()->registerObject(EventModel::Recon::RecMdcHitCol, hitList);
00436   if(!sc.isSuccess()) {
00437     log << MSG::FATAL << " Could not register RecMdcHit collection" <<endreq;
00438   }
00439   int trackId = trackList->size();
00440 
00441   if(m_debug>1)std::cout<<__FILE__<<" "<<__LINE__<<" STORED"<<  std::endl;
00442   MdcTrack mdcTrack(aTrack);//aTrack have been deleted in ~MdcTrack() 
00443 
00444   //tkStat:  0,PatRec; 1,MdcxReco; 2,Tsf; 3,CurlFinder; -1,Combined cosmic
00445   int tkStat = -1;
00446   mdcTrack.storeTrack(trackId, trackList, hitList, tkStat);
00447 
00448   //store hits on track
00449   RecMdcTrackCol::iterator it = trackList->begin();
00450   HitRefVec hl = (*it)->getVecHits();
00451   HitRefVec hitRefVec;
00452   HitRefVec::iterator itHit = hl.begin();   
00453   for (;itHit!=hl.end();++itHit){
00454     RecMdcHit* h = *itHit;
00455     SmartRef<RecMdcHit> refHit(h);
00456     hitRefVec.push_back(refHit);
00457   }
00458 
00459   //store skipped hits
00460   HitRefVec::iterator itHitSkipped = skippedHits.begin();   
00461   for (;itHitSkipped!=skippedHits.end();++itHitSkipped){
00462     hitList->push_back(*itHitSkipped);
00463     SmartRef<RecMdcHit> refHitSkipped(*itHitSkipped);
00464     hitRefVec.push_back(refHitSkipped);
00465   }
00466   (*it)->setVecHits(hitRefVec);
00467 
00468   if(m_debug) std::cout<<__FILE__<<" "<<__LINE__<<" stored "<< (*it)->getVecHits().size()<<  std::endl;
00469 }


Public Member Functions
StatusCode	beginRun ()
StatusCode	beginRun ()
void	clearRecMdcTrackHit ()
void	clearRecMdcTrackHit ()
StatusCode	execute ()
StatusCode	execute ()
StatusCode	finalize ()
StatusCode	finalize ()
StatusCode	initialize ()
StatusCode	initialize ()
	MdcxCosmicSewer (const std::string &name, ISvcLocator *pSvcLocator)
	MdcxCosmicSewer (const std::string &name, ISvcLocator *pSvcLocator)
void	MdcxHitsToHots (TrkHitList *trks, HitRefVec &hits, HitRefVec &skiped)
void	MdcxHitsToHots (TrkHitList *trks, HitRefVec &hits, HitRefVec &skiped)
void	store (TrkRecoTrk *tk, HitRefVec &skip)
void	store (TrkRecoTrk *tk, HitRefVec &skip)
virtual	~MdcxCosmicSewer ()
virtual	~MdcxCosmicSewer ()
Private Attributes
BField *	m_bfield
BField *	m_bfield
double	m_bunchT0
TrkContextEv *	m_context
TrkContextEv *	m_context
std::vector< float >	m_cosmicSewPar
std::vector< float >	m_cosmicSewPar
bool	m_cosmicSewSkip
bool	m_countPropTime
NTuple::Item< double >	m_csmcD0
NTuple::Item< double >	m_csmcD0
NTuple::Item< double >	m_csmcOmega
NTuple::Item< double >	m_csmcOmega
NTuple::Item< double >	m_csmcPhi0
NTuple::Item< double >	m_csmcPhi0
NTuple::Item< double >	m_csmcPt
NTuple::Item< double >	m_csmcPt
NTuple::Item< double >	m_csmcTanl
NTuple::Item< double >	m_csmcTanl
NTuple::Item< double >	m_csmcZ0
NTuple::Item< double >	m_csmcZ0
bool	m_debug
bool	m_doSag
const MdcDetector *	m_gm
const MdcDetector *	m_gm
bool	m_hist
const MdcCalibFunSvc *	m_mdcCalibFunSvc
const MdcCalibFunSvc *	m_mdcCalibFunSvc
MdcCheckUtil *	m_mdcCheckUtil
MdcCheckUtil *	m_mdcCheckUtil
IMagneticFieldSvc *	m_pIMF
IMagneticFieldSvc *	m_pIMF
RawDataProviderSvc *	m_rawDataProviderSvc
RawDataProviderSvc *	m_rawDataProviderSvc
NTuple::Tuple *	m_xtupleCsmcSew
NTuple::Tuple *	m_xtupleCsmcSew

MdcxCosmicSewer Class Reference

Public Member Functions

Private Attributes

Constructor & Destructor Documentation

Member Function Documentation

Member Data Documentation