/home/bes3soft/bes3soft/Boss/7.0.2/dist/7.0.2/Reconstruction/MdcPatRec/MdcxReco/MdcxReco-00-01-59/src/MdcxCosmicSewer.cxx

Go to the documentation of this file.
00001 //--------------------------------------------------------------------------
00002 // File and Version Information:
00003 //      $Id: MdcxCosmicSewer.cxx,v 1.11 2012/07/20 05:48:16 zhangy Exp $
00004 //
00005 // Description:
00006 //      Class MdcxCosmicSewer
00007 //
00008 // Environment:
00009 //      Software developed for the BaBar Detector at the SLAC B-Factory.
00010 //
00011 // Author List:
00012 //      Steve Wagner                    Original Author
00013 //      Zhang Yao(zhangyao@ihep.ac.cn)
00014 //
00015 // Copyright Information:
00016 //      Copyright (C) 1997              BEPCII
00017 // 
00018 // History:
00019 //      Migration for BESIII MDC
00020 //
00021 //------------------------------------------------------------------------
00022 
00023 //-----------------------
00024 // This Class's Header --
00025 //-----------------------
00026 #include "MdcxReco/MdcxCosmicSewer.h"
00027 
00028 //-------------------------------
00029 // Collaborating Class Headers --
00030 //-------------------------------
00031 #include "GaudiKernel/MsgStream.h"
00032 #include "GaudiKernel/AlgFactory.h"
00033 #include "EventModel/EventHeader.h"
00034 #include "EvTimeEvent/RecEsTime.h"
00035 #include "Identifier/MdcID.h"
00036 #include "Identifier/Identifier.h"
00037 #include "GaudiKernel/SmartDataPtr.h"
00038 #include "MdcCalibFunSvc/IMdcCalibFunSvc.h"
00039 #include "RawDataProviderSvc/IRawDataProviderSvc.h"
00040 #include "RawDataProviderSvc/MdcRawDataProvider.h"
00041 #include "MdcRawEvent/MdcDigi.h"
00042 #include "TrkBase/TrkFit.h"
00043 #include "TrkFitter/TrkHelixMaker.h"
00044 #include "TrkFitter/TrkLineMaker.h"
00045 #include "MdcxReco/Mdcxprobab.h"
00046 #include "MdcData/MdcHit.h"
00047 #include "MdcData/MdcRecoHitOnTrack.h"
00048 #include "MdcTrkRecon/MdcTrack.h"
00049 #include "TrkBase/TrkHitList.h"
00050 #include "MdcPrintSvc/IMdcPrintSvc.h"
00051 #include "TrkBase/TrkExchangePar.h"
00052 #include "MdcGeom/EntranceAngle.h"
00053 
00054 //----------------
00055 // Constructors --
00056 //----------------
00057 
00058 MdcxCosmicSewer::MdcxCosmicSewer(const std::string& name, ISvcLocator* pSvcLocator) :
00059   Algorithm(name, pSvcLocator) 
00060 {
00061 
00062   declareProperty("debug",          m_debug = false);
00063   declareProperty("hist",           m_hist  = false);
00064 
00065   declareProperty("doSag",          m_doSag = false);
00066 
00067   declareProperty("cosmicSewPar",   m_cosmicSewPar);
00068   declareProperty("cosmicSewSkip",  m_cosmicSewSkip=false);
00069 
00070   declareProperty("countPropTime",  m_countPropTime = true);
00071   declareProperty("lineFit",        m_lineFit = false);
00072 }
00073 
00074 //--------------
00075 // Destructor --
00076 //--------------
00077 MdcxCosmicSewer::~MdcxCosmicSewer() {
00078   delete m_bfield;
00079   delete m_context;
00080 }
00081 
00082 StatusCode MdcxCosmicSewer::beginRun(){  
00083   //Initailize MdcDetector
00084   m_gm = MdcDetector::instance(m_doSag); 
00085   if(NULL == m_gm) return StatusCode::FAILURE;
00086 
00087   return StatusCode::SUCCESS;
00088 }
00089 
00090 //--------------
00091 // Operations --
00092 //--------------
00093 StatusCode MdcxCosmicSewer::initialize(){  
00094   MsgStream log(msgSvc(), name());
00095   StatusCode sc;
00096   log << MSG::INFO << "in initialize()" << endreq;      
00097 
00098   i_evt=0;
00099   m_nSewed = 0;
00100 
00101   //Initailize magnetic filed 
00102   sc = service ("MagneticFieldSvc",m_pIMF);
00103   if(sc!=StatusCode::SUCCESS) {
00104     log << MSG::ERROR << "Unable to open Magnetic field service"<<endreq;
00105   }
00106   m_bfield = new BField(m_pIMF);
00107   log << MSG::INFO << "field z = "<<m_bfield->bFieldNominal()<< endreq;
00108   m_context = new TrkContextEv(m_bfield);
00109 
00110 
00111   // Get MdcCalibFunSvc
00112   IMdcCalibFunSvc* imdcCalibSvc; 
00113   sc = service ("MdcCalibFunSvc", imdcCalibSvc);
00114   m_mdcCalibFunSvc = dynamic_cast<MdcCalibFunSvc*>(imdcCalibSvc);
00115   if ( sc.isFailure() ){
00116     log << MSG::FATAL << "Could not load MdcCalibFunSvc!" << endreq;
00117     return StatusCode::FAILURE;
00118   }
00119 
00120   // Get RawDataProviderSvc
00121   IRawDataProviderSvc* iRawDataProvider; 
00122   sc = service ("RawDataProviderSvc", iRawDataProvider);
00123   if ( sc.isFailure() ){
00124     log << MSG::FATAL << "Could not load RawDataProviderSvc!" << endreq;
00125     return StatusCode::FAILURE;
00126   }
00127   m_rawDataProviderSvc = dynamic_cast<RawDataProviderSvc*>(iRawDataProvider);
00128 
00129   //Initailize MdcUtilitySvc
00130   IMdcUtilitySvc * imdcUtilitySvc;
00131   sc = service ("MdcUtilitySvc", imdcUtilitySvc);
00132   m_mdcUtilitySvc = dynamic_cast<MdcUtilitySvc*> (imdcUtilitySvc);
00133   if ( sc.isFailure() ){
00134     log << MSG::FATAL << "Could not load MdcUtilitySvc!" << endreq;
00135     return StatusCode::FAILURE;
00136   }
00137 
00138   //Initailize MdcPrintSvc
00139   IMdcPrintSvc* imdcPrintSvc;
00140   sc = service ("MdcPrintSvc", imdcPrintSvc);
00141   m_mdcPrintSvc = dynamic_cast<MdcPrintSvc*> (imdcPrintSvc);
00142   if ( sc.isFailure() ){
00143     log << MSG::FATAL << "Could not load MdcPrintSvc!" << endreq;
00144     return StatusCode::FAILURE;
00145   }
00146 
00147   if(m_hist){
00148     //book ntuple
00149     NTuplePtr ntCsmcSew(ntupleSvc(), "MdcxReco/csmc");
00150     if ( ntCsmcSew ) { m_xtupleCsmcSew = ntCsmcSew;}
00151     else {
00152       m_xtupleCsmcSew = ntupleSvc()->book ("MdcxReco/csmc", CLID_ColumnWiseTuple, "MdcxReco reconsturction results");
00153       if ( m_xtupleCsmcSew ) {
00154         sc = m_xtupleCsmcSew->addItem ("dD0",      m_csmcD0);
00155         sc = m_xtupleCsmcSew->addItem ("dPhi0",    m_csmcPhi0);
00156         sc = m_xtupleCsmcSew->addItem ("dZ0",      m_csmcZ0);
00157         sc = m_xtupleCsmcSew->addItem ("dOmega",   m_csmcOmega);
00158         sc = m_xtupleCsmcSew->addItem ("dPt",      m_csmcPt);
00159         sc = m_xtupleCsmcSew->addItem ("dTanl",    m_csmcTanl);
00160       }else{
00161         log << MSG::FATAL << "Could not book MdcxReco/csmc!" << endreq;
00162         return StatusCode::FAILURE;
00163       }
00164     }//end of book
00165   }
00166 
00167   return StatusCode::SUCCESS;
00168 }
00169 
00170 StatusCode MdcxCosmicSewer::execute() {
00171   MsgStream log(msgSvc(), name());
00172   log << MSG::INFO << "in execute()" << endreq; 
00173 
00174 
00175   setFilterPassed(false);
00176   // Get event No.
00177   SmartDataPtr<Event::EventHeader> evtHead(eventSvc(),"/Event/EventHeader");
00178   if (!evtHead) {
00179     log << MSG::FATAL<< "Could not retrieve event header" << endreq;
00180     return StatusCode::FAILURE;
00181   } 
00182   long t_runNo = evtHead->runNumber();
00183   long t_evtNo = evtHead->eventNumber();
00184   if (m_debug) std::cout << "sew "<<++i_evt<<" run "<<t_runNo<<" evt " << t_evtNo << std::endl;
00185 
00186 
00187   // Get bunch time t0 (ns)
00188   m_bunchT0 = -999.;
00189   SmartDataPtr<RecEsTimeCol> aevtimeCol(eventSvc(),"/Event/Recon/RecEsTimeCol");
00190   if (!aevtimeCol || aevtimeCol->size()==0) {
00191     log << MSG::WARNING<< "Could not find RecEsTimeCol"<< endreq;
00192     return StatusCode::SUCCESS;
00193   }
00194   RecEsTimeCol::iterator iter_evt = aevtimeCol->begin();
00195   for(; iter_evt!=aevtimeCol->end(); iter_evt++){
00196     m_bunchT0 =  (*iter_evt)->getTest();
00197   }
00198 
00199   //read event data 
00200   SmartDataPtr<RecMdcTrackCol> recMdcTrackCol(eventSvc(),EventModel::Recon::RecMdcTrackCol);
00201   SmartDataPtr<RecMdcHitCol> recMdcHitCol(eventSvc(),EventModel::Recon::RecMdcHitCol);
00202   if (!recMdcTrackCol || ! recMdcHitCol) return StatusCode::SUCCESS; 
00203   if (2!=recMdcTrackCol->size()){
00204     clearRecMdcTrackHit();
00205     return StatusCode::SUCCESS;
00206   }
00207   if(m_debug)std::cout<<name()<<" nTk==2 begin sewing"<<  std::endl;
00208 
00209 
00210   //get best track
00211   double dParam[5]={0.,0.,0.,0.,0.};
00212   float bprob =0.;
00213   float chisum =0.; 
00214   RecMdcTrackCol::iterator besthel;
00215 
00216   int besthelId = -999;
00217   RecMdcTrackCol::iterator it = recMdcTrackCol->begin();
00218   for (int iTk=0; it != recMdcTrackCol->end(); it++,iTk++) {
00219     RecMdcTrack* tk = *it;
00220     double Bz = m_bfield->bFieldNominal();
00221     double d0 = -tk->helix(0);    //cm
00222     double phi0 = tk->helix(1)+ CLHEP::halfpi;
00223     double omega = Bz*tk->helix(2)/-333.567;
00224     double z0 = tk->helix(3);    //cm
00225     double tanl = tk->helix(4);    //cm
00226 
00227     if(phi0>Constants::pi)phi0-=Constants::twoPi;
00228     if(phi0<-Constants::pi)phi0+=Constants::twoPi;
00229 
00230     if(m_debug)std::cout<<__FILE__<<" "<<__LINE__<< "tk"<<iTk
00231       <<"("<<d0<<","<<phi0<<"," <<","<<omega<<","<<z0<<","<<tanl<<")"<<std::endl;
00232 
00233     if(iTk==0){//store param1
00234       dParam[0]=d0;
00235       dParam[1]=phi0;
00236       dParam[2]=omega;
00237       dParam[3]=z0;
00238       dParam[4]=tanl;
00239     }else{//calc. diff between param1 and param2
00240       dParam[0] += d0;
00241       dParam[1] -= (phi0+Constants::pi);
00242       dParam[2] -= omega;
00243       dParam[3] -= z0;
00244       dParam[4] += tanl;
00245     }
00246 
00247     if(phi0<0) {
00248       besthelId = iTk;
00249       besthel = it;
00250     }
00251     float bchisq = tk->chi2();
00252     int   bndof  = tk->ndof();
00253     bprob=Mdcxprobab(bndof,bchisq);
00254     chisum += bchisq;
00255     //if (bprob > bestprob){
00256     //bestprob = bprob; 
00257     //besthel = it;
00258     //}
00259   }//end for recMdcTrackCol
00260 
00261   if(besthelId == -999) return StatusCode::SUCCESS;
00262 
00263 
00264   TrkHelixMaker helixfactory;
00265   TrkLineMaker linefactory; 
00266 
00267   if(m_debug){
00268     std::cout<<__FILE__<<" param diff: " << "\n D0 " << dParam[0]
00269       << "\n Phi0 " << dParam[1] << "\n Omega " << dParam[2]
00270       << "\n Z0 " << dParam[3] << "\n Tanl " << dParam[4] << std::endl;
00271   }
00272 
00273   if(m_hist){
00274     m_csmcD0= dParam[0];
00275     m_csmcPhi0=dParam[1];
00276     m_csmcOmega=dParam[2];
00277     m_csmcZ0=dParam[3];
00278     m_csmcTanl=dParam[4];
00279     m_xtupleCsmcSew->write();
00280   }
00281 
00282   //get track 5 parameters
00283   if( (fabs(dParam[0])>m_cosmicSewPar[0]) ||
00284       (fabs(dParam[1])>m_cosmicSewPar[1]) || 
00285       (fabs(dParam[2])>m_cosmicSewPar[2]) ||
00286       (fabs(dParam[3])>m_cosmicSewPar[3]) ||
00287       (fabs(dParam[4])>m_cosmicSewPar[4]) ){
00288     if(m_debug)std::cout<<__FILE__<<" 2 trk param not satisfied skip!"<<std::endl;
00289     if(m_debug){
00290       if(fabs(dParam[0])>m_cosmicSewPar[0]) std::cout<<" cut by d0 "<<std::endl;
00291       if(fabs(dParam[1])>m_cosmicSewPar[1]) std::cout<<" cut by phi0 "<<std::endl;
00292       if(fabs(dParam[2])>m_cosmicSewPar[2]) std::cout<<" cut by omega "<<std::endl;
00293       if(fabs(dParam[3])>m_cosmicSewPar[3]) std::cout<<" cut by z0"<<std::endl;
00294       if(fabs(dParam[4])>m_cosmicSewPar[4]) std::cout<<" cut by tanl"<<std::endl;
00295     }
00296     clearRecMdcTrackHit();
00297     return StatusCode::SUCCESS;
00298   }
00299 
00300   //got 2 sew-able trks, get helix param
00301   RecMdcTrack* tk = *besthel;
00302   double Bz = m_bfield->bFieldNominal();
00303   double d0 = -tk->helix(0);    //cm
00304   double phi0 = tk->helix(1)+ CLHEP::halfpi;
00305   double omega = Bz*tk->helix(2)/-333.567;
00306   double z0 = tk->helix(3);    //cm
00307   double tanl = tk->helix(4);    //cm
00308 
00309   TrkExchangePar tt(d0,phi0,omega,z0,tanl);
00310   if(m_debug)std::cout<<__FILE__<<" best helix: No."<<tk->trackId()<<" Pat param=("<<d0<<" "<<phi0
00311     <<" "<<omega<<" "<<z0<<" "<<tanl<<")"<<  std::endl;
00312 
00313 
00314   //new track
00315   TrkRecoTrk* newTrack;
00316   if(m_lineFit){
00317     newTrack = linefactory.makeTrack(tt, chisum, *m_context, m_bunchT0*1.e-9);
00318     linefactory.setFlipAndDrop(*newTrack, true, true);
00319   }else{
00320     newTrack = helixfactory.makeTrack(tt, chisum, *m_context, m_bunchT0*1.e-9);
00321     helixfactory.setFlipAndDrop(*newTrack, true, true);
00322   }
00323 
00324 
00325   // combine hit list
00326   HitRefVec skippedHits; 
00327   it = recMdcTrackCol->begin();
00328   for (int iTk=0; it != recMdcTrackCol->end(); it++,iTk++) {
00329     RecMdcTrack* tk = *it;
00330     HitRefVec hits = tk->getVecHits();
00331     HepVector helixPatPar = m_mdcUtilitySvc->besPar2PatPar(tk->helix());
00332     MdcxHitsToHots(helixPatPar,newTrack,hits,skippedHits);
00333   }
00334 
00335 
00336   //------------------------------------
00337   //do fit
00338   //------------------------------------
00339   TrkErrCode err=newTrack->hits()->fit();
00340   const TrkFit* newFit = newTrack->fitResult(); 
00341 
00342 
00343   if (!newFit) {
00344     // Fit bad
00345     if(m_debug)std::cout<<__FILE__<<" sew fit failed!!!"<<  std::endl;
00346     HitRefVec::iterator it= skippedHits.begin();   
00347     for (;it!=skippedHits.end();++it){ delete it->data();}
00348     delete newTrack;
00349     clearRecMdcTrackHit();
00350   } else {
00351     //------------------------------------
00352     // Fit good
00353     //------------------------------------
00354     if(m_lineFit){
00355       linefactory.setFlipAndDrop(*newTrack, false, false);
00356     }else{
00357       helixfactory.setFlipAndDrop(*newTrack, false, false);
00358     }
00359     if(m_debug){
00360       std::cout<<__FILE__<<" sew cosmic fit good "<<  std::endl;
00361       newTrack->print(std::cout);
00362     }
00363 
00364     SmartDataPtr<MdcHitCol> m_hitCol(eventSvc(), "/Event/Hit/MdcHitCol");
00365     if (!m_hitCol){
00366       m_hitCol= new MdcHitCol;
00367       StatusCode sc = eventSvc()->registerObject("/Event/Hit/MdcHitCol",m_hitCol);
00368       if(!sc.isSuccess()) {
00369         std::cout<< " Could not register MdcHitCol" <<endreq;
00370         return StatusCode::FAILURE;
00371       }
00372     }
00373     uint32_t getDigiFlag = 0;
00374     MdcDigiVec mdcDigiVec = m_rawDataProviderSvc->getMdcDigiVec(getDigiFlag);
00375     const MdcDigi* m_digiMap[43][288];
00376     MdcDigiVec::iterator iter = mdcDigiVec.begin();
00377     for (; iter != mdcDigiVec.end(); iter++ ) {
00378       const MdcDigi* aDigi = *iter;
00379       const Identifier id= aDigi->identify();
00380       int layer = MdcID::layer(id);
00381       int wire = MdcID::wire(id);
00382       m_digiMap[layer][wire] = aDigi;
00383     }
00384 
00385     //calc. doca of skipped hits
00386     HitRefVec::iterator itHitSkipped = skippedHits.begin();   
00387     for (;itHitSkipped!=skippedHits.end();++itHitSkipped){
00388       RecMdcHit* h = *itHitSkipped;
00389       Identifier id(h->getMdcId());
00390       int layer = MdcID::layer(id);
00391       int wire  = MdcID::wire(id);
00392       double fltLen = 0.;
00393       HepVector helix = newFit->helix(fltLen).params();
00394       HepSymMatrix err = newFit->helix(fltLen).covariance();
00395       double docaFltLen = h->getFltLen();
00396       //double doca = m_mdcUtilitySvc->docaPatPar(layer,wire,docaFltLen,helix,err);//yzhang 2012-07-19 
00397       double doca = m_mdcUtilitySvc->docaPatPar(layer,wire,helix,err);
00398       double newDoca = fabs(doca);
00399       int flagLR = h->getFlagLR();
00400       if ( flagLR == 0 ){ newDoca = -fabs(doca); }//left driftDist <0
00401       h->setDoca(newDoca);
00402       if(m_debug>=3)std::cout<<"("<<layer<<","<<wire<<") new doca "<<newDoca<<" resid "<<fabs(fabs(newDoca)-fabs(h->getDriftDistLeft()))<<std::endl;
00403       if(m_debug>=3&&fabs(fabs(newDoca)-fabs(h->getDriftDistLeft()))>0.02) std::cout<<" bad resid  "<<fabs(fabs(newDoca)-fabs(h->getDriftDistLeft()))<<" doca  "<<newDoca<<" dd  "<<h->getDriftDistLeft()<<std::endl;
00404 
00405       //calc drift dist of skipped hit
00406       MdcHit *thehit = new MdcHit(m_digiMap[layer][wire], m_gm);
00407       thehit->setCosmicFit(true);
00408       thehit->setCalibSvc(m_mdcCalibFunSvc);
00409       thehit->setCountPropTime(m_countPropTime);
00410       m_hitCol->push_back(thehit);
00411 
00412       HepPoint3D poca; Hep3Vector tempDir;
00413       getInfo(helix,0,poca,tempDir);
00414       if(m_debug>3)std::cout<<"track poca  "<<poca<<std::endl;
00415       HepPoint3D pos; Hep3Vector dir;
00416       getInfo(helix,docaFltLen,pos,dir);
00417       //std::cout<<" pos  "<<pos<<" dir  "<<dir<<"   "<<std::endl;
00418       if(pos.y()>poca.y()){
00419         int wireAmb = patAmbig(flagLR);
00420         //double tofold = (docaFltLen/Constants::c + m_bunchT0)*1.e-9;
00421         double tof = docaFltLen/Constants::c + (m_bunchT0)*1.e-9;
00422         //std::cout<<layer<<" "<<wire<<" tofold  "<<tofold*1.e9-m_bunchT0 <<" tof  "<<tof*1.e9-m_bunchT0<<" diff  "<<(tof-tofold)*1.e9<<std::endl;
00423         //tof = tofold;
00424         double eAngle = EntranceAngle(dir.phi() - thehit->phi(pos.z()));
00425         double dAngle = Constants::pi/2 - 1.*dir.theta();
00426         double z = pos.z();
00427         double dist = thehit->driftDist(tof, wireAmb, eAngle, dAngle, z);
00428         double sigma = thehit->sigma(dist,wireAmb,eAngle,dAngle,z);
00429 
00430         if(m_debug>3&&fabs(fabs(h->getDoca())-fabs(dist))>0.02)
00431           std::cout<<"("<<layer<<","<<wire<<") old doca "<<h->getDoca()<<" old ddr "<<h->getDriftDistRight()
00432             <<" new dd "<<dist <<" newAmbig "<<wireAmb
00433             //<<" tof "<<tof <<" z "<<z <<" eAngle "<<eAngle <<" dTime "<<thehit->driftTime(tof,z)
00434             <<" old sigma  "<<h->getErrDriftDistLeft()
00435             <<" new sigma  "<<sigma<<"   "<<std::endl;
00436         h->setDriftDistLeft(-1.*abs(dist));
00437         h->setDriftDistRight(abs(dist));
00438         h->setErrDriftDistLeft(sigma);
00439         h->setErrDriftDistRight(sigma);
00440       }
00441     }
00442 
00443     //-------------------------------------------------------
00444     //Store TDS
00445     //-------------------------------------------------------
00446     store(newTrack, skippedHits);//newTrack have been deleted
00447 
00448     setFilterPassed(true);
00449     m_nSewed++;
00450   }// end if newFit
00451 
00452   if(m_debug) { m_mdcPrintSvc->printRecMdcTrackCol(); } 
00453 
00454   return StatusCode::SUCCESS;
00455 }
00456 
00457 StatusCode MdcxCosmicSewer::finalize() {
00458   MsgStream log(msgSvc(), name());
00459   log << MSG::INFO << "in finalize()" << endreq;
00460   std::cout<<name()<<" after sewed got "<<m_nSewed<<" cosmic tracks"<<std::endl;
00461   return StatusCode::SUCCESS;
00462 }
00463 
00464 void MdcxCosmicSewer::MdcxHitsToHots(HepVector& helix, TrkRecoTrk* newTrack, HitRefVec& recMdcHits, HitRefVec& skippedHits) {
00465 
00466   TrkHitList* trkHitList = newTrack->hits();
00467   //store new MdcHit
00468   SmartDataPtr<MdcHitCol> m_hitCol(eventSvc(), "/Event/Hit/MdcHitCol");
00469   if (!m_hitCol){
00470     m_hitCol= new MdcHitCol;
00471     StatusCode sc = eventSvc()->registerObject("/Event/Hit/MdcHitCol",m_hitCol);
00472     if(!sc.isSuccess()) {
00473       std::cout<< " Could not register MdcHitCol" <<endreq;
00474       return;
00475     }
00476   } 
00477 
00478   //get MdcDigi pointer
00479   uint32_t getDigiFlag = 0;
00480   MdcDigiVec mdcDigiVec = m_rawDataProviderSvc->getMdcDigiVec(getDigiFlag); 
00481   if (0 == mdcDigiVec.size()){
00482     std::cout << " No hits in MdcDigiVec" << std::endl;
00483     return; 
00484   }
00485   const MdcDigi* m_digiMap[43][288];
00486   MdcDigiVec::iterator iter = mdcDigiVec.begin();
00487   for (; iter != mdcDigiVec.end(); iter++ ) {
00488     const MdcDigi* aDigi = *iter;
00489     const Identifier id= aDigi->identify();
00490     int layer = MdcID::layer(id);
00491     int wire = MdcID::wire(id);
00492     m_digiMap[layer][wire] = aDigi;
00493   }
00494 
00495 
00496   HepPoint3D poca; Hep3Vector tempDir;
00497   getInfo(helix,0,poca,tempDir);
00498   if(m_debug>3)std::cout<<"track poca  "<<poca<<std::endl;
00499   // new MdcRecoHitOnTracks
00500   HitRefVec::iterator itHit = recMdcHits.begin();   
00501   for (;itHit!=recMdcHits.end();++itHit){
00502     RecMdcHit* h = *itHit;
00503     int layer = MdcID::layer(h->getMdcId());
00504     int wire = MdcID::wire(h->getMdcId());
00505     // new fltLen and ambig
00506     double newFltLen = h->getFltLen();
00507     HepPoint3D pos; Hep3Vector dir;
00508     getInfo(helix,h->getFltLen(),pos,dir);
00509     int newAmbig = patAmbig(h->getFlagLR());
00510     if(pos.y()>poca.y()){//yzhang TEMP 2012-07-17 
00511       newFltLen *= -1.;
00512       newAmbig *= -1;
00513       if(m_debug>3) std::cout<<" Up track, change sign of fltLen  "<<std::endl;
00514     }
00515     int newFlagLR = bes3FlagLR(newAmbig);
00516     // calc. position of this point
00517     if(m_cosmicSewSkip && layer<20){
00518       RecMdcHit* newSkippedHit = new RecMdcHit();
00519       newSkippedHit->setDriftDistLeft(h->getDriftDistLeft());
00520       newSkippedHit->setDriftDistRight(h->getDriftDistRight());
00521       newSkippedHit->setErrDriftDistLeft(h->getErrDriftDistLeft());
00522       newSkippedHit->setErrDriftDistRight(h->getErrDriftDistRight());
00523       newSkippedHit->setChisqAdd(h->getChisqAdd());
00524       newSkippedHit->setFlagLR(newFlagLR);
00525       newSkippedHit->setStat(h->getStat());
00526       newSkippedHit->setMdcId(h->getMdcId());
00527       newSkippedHit->setTdc(h->getTdc());
00528       newSkippedHit->setAdc(h->getAdc());
00529       newSkippedHit->setDriftT(h->getDriftT());
00530       newSkippedHit->setDoca(999.);
00531       //newSkippedHit->setDoca(h->getDoca());
00532       newSkippedHit->setEntra(h->getEntra());
00533       newSkippedHit->setZhit(h->getZhit());
00534       newSkippedHit->setFltLen(newFltLen);
00535       skippedHits.push_back(newSkippedHit);
00536     }else{
00537       MdcHit *thehit = new MdcHit(m_digiMap[layer][wire], m_gm); 
00538       thehit->setCosmicFit(true);
00539       thehit->setCalibSvc(m_mdcCalibFunSvc);
00540       thehit->setCountPropTime(m_countPropTime);
00541       m_hitCol->push_back(thehit);
00542 
00543       MdcRecoHitOnTrack temp(*thehit, newAmbig, m_bunchT0);//m_bunchT0 nano second here
00544       MdcHitOnTrack* newhot = &temp;
00545       newhot->setFltLen(newFltLen);
00546       if(m_debug>3) std::cout<<name()<<" ("<<layer<<","<<wire<<") fltLen "<<h->getFltLen()<<" newFltLen "<<newFltLen<<" newAmbig "<<newAmbig<<" pos.y() "<<pos.y()<<std::endl;
00547 
00548       trkHitList->appendHot(newhot);
00549     }
00550   }//end of loop hits
00551 }
00552 
00553 void MdcxCosmicSewer::store(TrkRecoTrk* aTrack, HitRefVec& skippedHits) {
00554   StatusCode sc;
00555   MsgStream log(msgSvc(), name());
00556 
00557   assert (aTrack != NULL);
00558 
00559   //IDataManagerSvc *dataManSvc;
00560   SmartIF<IDataManagerSvc> dataManSvc(eventSvc()); 
00561   DataObject *aTrackCol;
00562   eventSvc()->findObject("/Event/Recon/RecMdcTrackCol",aTrackCol);
00563   if(aTrackCol != NULL) {
00564     dataManSvc->clearSubTree("/Event/Recon/RecMdcTrackCol");
00565     eventSvc()->unregisterObject("/Event/Recon/RecMdcTrackCol");
00566   }
00567 
00568   DataObject *aRecHitCol;
00569   eventSvc()->findObject("/Event/Recon/RecMdcHitCol",aRecHitCol);
00570   if(aRecHitCol != NULL) {
00571     dataManSvc->clearSubTree("/Event/Recon/RecMdcHitCol");
00572     eventSvc()->unregisterObject("/Event/Recon/RecMdcHitCol");
00573   }
00574 
00575   SmartDataPtr<RecMdcTrackCol> trackList(eventSvc(),EventModel::Recon::RecMdcTrackCol);
00576   if (!trackList) {
00577     RecMdcTrackCol *trackList = new RecMdcTrackCol;
00578     sc =  eventSvc()->registerObject(EventModel::Recon::RecMdcTrackCol, trackList);
00579     if(!sc.isSuccess()) {
00580       log << MSG::FATAL << " Could not register RecMdcTrack collection" <<endreq;
00581     }
00582   }
00583 
00584   SmartDataPtr<RecMdcHitCol> hitList(eventSvc(),EventModel::Recon::RecMdcHitCol);
00585   if (!hitList) {
00586     hitList = new RecMdcHitCol;
00587     sc =  eventSvc()->registerObject(EventModel::Recon::RecMdcHitCol, hitList);
00588     if(!sc.isSuccess()) {
00589       log << MSG::FATAL << " Could not register RecMdcHit collection" <<endreq;
00590     }
00591   }
00592   int trackId = trackList->size() - 1;
00593 
00594   if(m_debug)std::cout<<__FILE__<<" "<<__LINE__<<" STORED"<<  std::endl;
00595   MdcTrack mdcTrack(aTrack);//aTrack have been deleted in ~MdcTrack() 
00596 
00597   //tkStat:  0,PatRec; 1,MdcxReco; 2,Tsf; 3,CurlFinder; -1,Combined cosmic
00598   int tkStat = -1;
00599   mdcTrack.storeTrack(trackId, trackList, hitList, tkStat);
00600 
00601   //store hits on track
00602   RecMdcTrackCol::iterator it = trackList->begin();
00603   HitRefVec hl = (*it)->getVecHits();
00604   HitRefVec hitRefVec;
00605   HitRefVec::iterator itHit = hl.begin();   
00606   int ihl=0;
00607   for (;itHit!=hl.end();++itHit){
00608     ihl++;
00609     RecMdcHit* h = *itHit;
00610     SmartRef<RecMdcHit> refHit(h);
00611     hitRefVec.push_back(refHit);
00612   }
00613 
00614   if(m_debug) std::cout<<__FILE__<<" "<<__LINE__<<"refHit stored "<< ihl<<  std::endl;
00615   //store skipped hits
00616   HitRefVec::iterator itHitSkipped = skippedHits.begin();   
00617   int iskipped=0;
00618   for (;itHitSkipped!=skippedHits.end();++itHitSkipped){
00619     iskipped++;
00620     (*itHitSkipped)->setTrkId(trackId);
00621     hitList->push_back(*itHitSkipped);
00622     SmartRef<RecMdcHit> refHitSkipped(*itHitSkipped);
00623     hitRefVec.push_back(refHitSkipped);
00624   }
00625   if(m_debug) std::cout<<__FILE__<<" "<<__LINE__<<"skippedHits stored "<< iskipped<<  std::endl;
00626   (*it)->setVecHits(hitRefVec);
00627 
00628   (*it)->setNhits((*it)->getVecHits().size());
00629   (*it)->setNster(-1);
00630   //omega==0, line fit :ndof = nhit-4
00631   //omega>0, helix fit :ndof = nhit-5
00632   if((*it)->helix(2)<0.00001)(*it)->setNdof((*it)->getNhits()-4);
00633   else (*it)->setNdof((*it)->getNhits()-5);
00634 
00635   if(m_debug) std::cout<<__FILE__<<" "<<__LINE__<<" stored "<< (*it)->getVecHits().size()<<  std::endl;
00636 
00637 }
00638 
00639 void MdcxCosmicSewer::clearRecMdcTrackHit(){
00640   MsgStream log(msgSvc(), name());
00641   StatusCode sc;
00642   //-------------------------------------------------------
00643   //clear TDS
00644   //-------------------------------------------------------
00645   //IDataManagerSvc *dataManSvc;
00646   SmartIF<IDataManagerSvc> dataManSvc(eventSvc()); 
00647   DataObject *aTrackCol;
00648   eventSvc()->findObject("/Event/Recon/RecMdcTrackCol",aTrackCol);
00649   if(aTrackCol != NULL) {
00650     dataManSvc->clearSubTree("/Event/Recon/RecMdcTrackCol");
00651     eventSvc()->unregisterObject("/Event/Recon/RecMdcTrackCol");
00652   }
00653   DataObject *aRecHitCol;
00654   eventSvc()->findObject("/Event/Recon/RecMdcHitCol",aRecHitCol);
00655   if(aRecHitCol != NULL) {
00656     dataManSvc->clearSubTree("/Event/Recon/RecMdcHitCol");
00657     eventSvc()->unregisterObject("/Event/Recon/RecMdcHitCol");
00658   }
00659   SmartDataPtr<RecMdcTrackCol> trackListTemp(eventSvc(),EventModel::Recon::RecMdcTrackCol);
00660   if (!trackListTemp) {
00661     RecMdcTrackCol *trackListTemp= new RecMdcTrackCol;
00662     sc =  eventSvc()->registerObject(EventModel::Recon::RecMdcTrackCol, trackListTemp);
00663     if(!sc.isSuccess()) {
00664       log << MSG::FATAL << " Could not register RecMdcTrack collection" <<endreq;
00665     }
00666   }
00667   SmartDataPtr<RecMdcHitCol> hitListTemp(eventSvc(),EventModel::Recon::RecMdcHitCol);
00668   if (!hitListTemp) {
00669     RecMdcHitCol *hitListTemp= new RecMdcHitCol;
00670     sc =  eventSvc()->registerObject(EventModel::Recon::RecMdcHitCol, hitListTemp);
00671     if(!sc.isSuccess()) {
00672       log << MSG::FATAL << " Could not register RecMdcHit collection" <<endreq;
00673     }
00674   }
00675 }
00676 
00677 void MdcxCosmicSewer::getInfo(HepVector helix, double fltLen, HepPoint3D& pos, Hep3Vector & dir){
00678   if(m_lineFit) return;//line fit FIXME yzhang 2012-07-19 
00679   double d0     = helix[0];
00680   double phi0   = helix[1];
00681   double omega  = helix[2];
00682   double z0     = helix[3];
00683   double tanDip = helix[4];
00684   double cDip = 1./sqrt(1.+tanDip*tanDip); 
00685   double sDip = tanDip * cDip;
00686   double phi00 = phi0;  // Don't normalize
00687   double ang = phi00 + cDip*fltLen*omega; 
00688   double cang = cos(ang);
00689   double sang = sin(ang);
00690   double sphi0 = sin(phi00);
00691   double cphi0 = cos(phi00);
00692   HepPoint3D referencePoint(0.,0.,0.);//yzhang 2012-07-17 TEMP
00693   double xt = (sang - sphi0)/omega - d0*sphi0 + referencePoint.x();
00694   double yt = -(cang - cphi0)/omega + d0*cphi0 + referencePoint.y();
00695   double zt = z0 + fltLen*sDip + referencePoint.z();
00696   pos.setX(xt);
00697   pos.setY(yt);
00698   pos.setZ(zt);
00699 
00700   dir.setX(cang * cDip);
00701   dir.setY(sang * cDip);
00702   dir.setZ(sDip);
00703 }
00704 
00705 int MdcxCosmicSewer::patAmbig(int bes3FlagLR){
00706   int ambig = 0;
00707   if ( bes3FlagLR == 0 ){ 
00708     ambig = 1; //left driftDist <0
00709   }else if( bes3FlagLR == 1){ 
00710     ambig = -1;//right driftDist >0
00711   }else if( bes3FlagLR == 2){ 
00712     ambig = 0; //don't know 
00713   }
00714   return ambig;
00715 }
00716 
00717 int MdcxCosmicSewer::bes3FlagLR(int patAmbig){
00718   int flagLR = 2;
00719   if ( patAmbig == 1 ){ 
00720     flagLR = 0; //left driftDist <0
00721   }else if( patAmbig == -1){ 
00722     flagLR = 1; //right driftDist >0
00723   }else if( patAmbig == 0){ 
00724     flagLR = 2; //don't know 
00725   }
00726   return flagLR;
00727 }

Generated on Tue Nov 29 23:13:37 2016 for BOSS_7.0.2 by  doxygen 1.4.7