#include <MdcxTrackFinder.h>
|
||||||||||||
|
00107 : 00108 Algorithm(name, pSvcLocator), m_mdcCalibFunSvc(0) 00109 { 00110 //input 00111 declareProperty("pdtFile", m_pdtFile = "pdt.table"); 00112 //debug control 00113 declareProperty("debug", m_debug= 0); 00114 declareProperty("hist", m_hist= 0); 00115 declareProperty("mcHist", m_mcHist = false); 00116 //cuts and control 00117 declareProperty("cresol", m_cresol = 0.013); 00118 00119 declareProperty("getDigiFlag", m_getDigiFlag = 0); 00120 declareProperty("maxMdcDigi", m_maxMdcDigi= 0); 00121 declareProperty("keepBadTdc", m_keepBadTdc= 0); 00122 declareProperty("dropHot", m_dropHot= 0); 00123 declareProperty("keepUnmatch", m_keepUnmatch= 0); 00124 declareProperty("salvageTrk", m_salvageTrk = false); 00125 declareProperty("dropMultiHotInLayer",m_dropMultiHotInLayer = false); 00126 declareProperty("dropTrkPt", m_dropTrkPt = -999.); 00127 declareProperty("d0Cut", m_d0Cut = 999.); 00128 declareProperty("z0Cut", m_z0Cut = 999.); 00129 00130 declareProperty("minMdcDigi", m_minMdcDigi = 0); 00131 declareProperty("countPropTime",m_countPropTime = true); 00132 declareProperty("addHitCut", m_addHitCut = 5.); 00133 declareProperty("dropHitsSigma",m_dropHitsSigma); 00134 declareProperty("helixFitCut", m_helixFitCut); 00135 declareProperty("minTrkProb", m_minTrkProb = 0.01); 00136 declareProperty("csmax4", m_csmax4 = 50.); 00137 declareProperty("csmax3", m_csmax3 = 1.); 00138 declareProperty("helixFitSigma",m_helixFitSigma= 5.); 00139 declareProperty("maxRcsInAddSeg",m_maxRcsInAddSeg= 50.); 00140 declareProperty("nSigAddHitTrk", m_nSigAddHitTrk= 5.); 00141 declareProperty("maxProca", m_maxProca= 0.6); 00142 declareProperty("doSag", m_doSag= false); 00143 declareProperty("lineFit", m_lineFit = false); 00144 declareProperty("cosmicFit", m_cosmicFit= false); 00145 }
|
|
|
00150 {
00151 delete m_bfield;
00152 }
|
|
||||||||||||
|
|
|
|
|
|
|
|
|
|
00154 {
00155 // Get Mdc Detector Geometry
00156 m_gm = MdcDetector::instance(m_doSag);
00157 if(NULL == m_gm) return StatusCode::FAILURE;
00158 MdcxHit::setMdcDetector(m_gm);
00159
00160 return StatusCode::SUCCESS;
00161 }
|
|
|
|
|
|
00788 {
00789 MsgStream log(msgSvc(), name());
00790 StatusCode sc;
00791 if(!m_hist) return;
00792
00793 g_poison = histoSvc()->book( "poison", "poison",43,0,42,288,0,288 );
00794 g_csmax4 = histoSvc()->book( "csmax4", "csmax4",100,0,100 );
00795 g_csmax3 = histoSvc()->book( "csmax3", "csmax3",50000,0,500000 );
00796 g_omegag = histoSvc()->book( "omegag", "omegag",1000 ,0.,1.);
00797 g_dPhiAU = histoSvc()->book( "dPhiAU", "dPhiAU",1000 ,0.,3.5);
00798 g_dPhiAU_0 = histoSvc()->book( "dPhiAU_0", "dPhiAU_0",1000 ,0.,3.5);
00799 g_dPhiAU_1 = histoSvc()->book( "dPhiAU_1", "dPhiAU_1",1000 ,0.,3.5);
00800 g_dPhiAU_5 = histoSvc()->book( "dPhiAU_5", "dPhiAU_5",1000 ,0.,3.5);
00801 g_dPhiAU_7 = histoSvc()->book( "dPhiAU_7", "dPhiAU_7",1000 ,0.,3.5);
00802 g_dPhiAV = histoSvc()->book( "dPhiAV", "dPhiAV",1000 ,0.,3.5);
00803 g_addSegPhi = histoSvc()->book( "addSegPhi", "addSegPhi",1000 ,0.,3.5);
00804 g_dPhiAV_0 = histoSvc()->book( "dPhiAV_0", "dPhiAV_0",1000 ,0.,3.5);
00805 g_dPhiAV_1 = histoSvc()->book( "dPhiAV_1", "dPhiAV_1",1000 ,0.,3.5);
00806 g_dPhiAV_6 = histoSvc()->book( "dPhiAV_6", "dPhiAV_6",1000 ,0.,3.5);
00807 g_dPhiAV_8 = histoSvc()->book( "dPhiAV_8", "dPhiAV_8",1000 ,0.,3.5);
00808 g_trkllmk = histoSvc()->book( "trkllmk", "trkllmk",43,0.,42);
00809 g_trklgood = histoSvc()->book( "trklgood", "trklgood",43,0.,42);
00810 g_trklcircle = histoSvc()->book( "trklcircle", "trklcircle",43,0.,42);
00811 g_trklhelix= histoSvc()->book( "trklhelix", "trklhelix",43,0.,42);
00812 g_trkldrop1= histoSvc()->book( "trkldrop1", "trkldrop1",43,0.,42);
00813 g_trkldrop2= histoSvc()->book( "trkldrop2", "trkldrop2",43,0.,42);
00814 g_trklappend1= histoSvc()->book( "trklappend1", "trklappend1",43,0.,42);
00815 g_trklappend2= histoSvc()->book( "trklappend2", "trklappend2",43,0.,42);
00816 g_trklappend3= histoSvc()->book( "trklappend3", "trklappend3",43,0.,42);
00817 //g_fitOmega= histoSvc()->book( "fitOmega", "fitOmega",200,0.,100.);
00818 g_trklfirstProb= histoSvc()->book( "trklfirstProb", "trklfirstProb",200,0.,2.);
00819 g_trkltemp= histoSvc()->book( "trkltemp", "trkltemp",200,-3.5,3.5);
00820
00821 g_trklproca= histoSvc()->book( "trklproca", "trklproca",100,0.,5.);
00822 g_trkle= histoSvc()->book( "trkle", "trkle",200,0.,0.025);
00823 g_trkld= histoSvc()->book( "trkld", "trkld",200,-1.2,1.2);
00824 g_trklel= histoSvc()->book( "trklel", "trklel",200,0.,0.025,43,0,43);
00825 g_trkldl= histoSvc()->book( "trkldl", "trkldl",200,-1.2,1.2,43,0,43);
00826 g_trkldoca= histoSvc()->book( "trkldoca", "trkldoca",200,-1.2,1.2);
00827 g_trkllayer= histoSvc()->book( "trkllayer", "trkllayer",43,0,43);
00828 g_dropHitsSigma= histoSvc()->book( "dropHitsSigma", "dropHitsSigma",43,0,43,100,0,11);
00829 g_addHitCut= histoSvc()->book( "addHitCut", "addHitCut",100,0,10);
00830 g_addHitCut2d= histoSvc()->book( "addHitCut2d", "addHitCut2d",43,0,43,100,0,10);
00831
00832 NTuplePtr nt1(ntupleSvc(), "FILE121/rec");
00833 if ( nt1 ) { m_xtuple1 = nt1;}
00834 else {
00835 m_xtuple1 = ntupleSvc()->book ("FILE121/rec", CLID_ColumnWiseTuple, "MdcxReco reconsturction results");
00836 if ( m_xtuple1 ) {
00837 sc = m_xtuple1->addItem ("t0", m_xt0);
00838 sc = m_xtuple1->addItem ("timing", m_xtiming);
00839 sc = m_xtuple1->addItem ("t0Stat", m_xt0Stat);
00840 sc = m_xtuple1->addItem ("t0Truth", m_xt0Truth);
00841
00842 sc = m_xtuple1->addItem ("nSlay", m_xnSlay);
00843 sc = m_xtuple1->addItem ("p", m_xp);
00844 sc = m_xtuple1->addItem ("pt", m_xpt);
00845 sc = m_xtuple1->addItem ("pz", m_xpz);
00846 sc = m_xtuple1->addItem ("d0", m_xd0);
00847 sc = m_xtuple1->addItem ("phi0", m_xphi0);
00848 sc = m_xtuple1->addItem ("cpa", m_xcpa);
00849 sc = m_xtuple1->addItem ("z0", m_xz0);
00850 sc = m_xtuple1->addItem ("tanl", m_xtanl);
00851 sc = m_xtuple1->addItem ("q", m_xq);
00852 sc = m_xtuple1->addItem ("pocax", m_xpocax);
00853 sc = m_xtuple1->addItem ("pocay", m_xpocay);
00854 sc = m_xtuple1->addItem ("pocaz", m_xpocaz);
00855
00856 sc = m_xtuple1->addItem ("evtNo", m_xevtNo);
00857 sc = m_xtuple1->addItem ("nSt", m_xnSt);
00858 sc = m_xtuple1->addItem ("nAct", m_xnAct);
00859 sc = m_xtuple1->addItem ("nDof", m_xnDof);
00860 sc = m_xtuple1->addItem ("chi2", m_xchi2);
00861 sc = m_xtuple1->addItem ("tkId", m_xtkId);
00862 sc = m_xtuple1->addItem ("nHit", m_xnHit, 0, 7000);//# of hit/rec track
00863 sc = m_xtuple1->addIndexedItem ("resid", m_xnHit, m_xresid);
00864 sc = m_xtuple1->addIndexedItem ("sigma", m_xnHit, m_xsigma);
00865 sc = m_xtuple1->addIndexedItem ("driftD", m_xnHit, m_xdriftD);
00866 sc = m_xtuple1->addIndexedItem ("driftT", m_xnHit, m_xdriftT);
00867 sc = m_xtuple1->addIndexedItem ("doca", m_xnHit, m_xdoca);
00868 sc = m_xtuple1->addIndexedItem ("entra", m_xnHit, m_xentra);
00869 sc = m_xtuple1->addIndexedItem ("ambig", m_xnHit, m_xambig);
00870 sc = m_xtuple1->addIndexedItem ("fltLen", m_xnHit, m_xfltLen);
00871 sc = m_xtuple1->addIndexedItem ("tof", m_xnHit, m_xtof);
00872 sc = m_xtuple1->addIndexedItem ("act", m_xnHit, m_xact);
00873 sc = m_xtuple1->addIndexedItem ("tdc", m_xnHit, m_xtdc);
00874 sc = m_xtuple1->addIndexedItem ("adc", m_xnHit, m_xadc);
00875 sc = m_xtuple1->addIndexedItem ("layer", m_xnHit, m_xlayer);
00876 sc = m_xtuple1->addIndexedItem ("wire", m_xnHit, m_xwire);
00877 sc = m_xtuple1->addIndexedItem ("x", m_xnHit, m_xx);
00878 sc = m_xtuple1->addIndexedItem ("y", m_xnHit, m_xy);
00879 sc = m_xtuple1->addIndexedItem ("z", m_xnHit, m_xz);
00880 } else { // did not manage to book the N tuple....
00881 log << MSG::ERROR << " Cannot book N-tuple:" << long(m_xtuple1) << endmsg;
00882 return;
00883 }
00884 }//end book of nt1
00885 NTuplePtr nt4(ntupleSvc(), "FILE121/evt");
00886 if ( nt4 ) { m_xtupleEvt = nt4;}
00887 else {
00888 m_xtupleEvt = ntupleSvc()->book ("FILE121/evt", CLID_ColumnWiseTuple, "MdcxReco event data");
00889 if ( m_xtupleEvt ) {
00890 sc = m_xtupleEvt->addItem ("evtNo", m_xt4EvtNo );
00891 sc = m_xtupleEvt->addItem ("nRecTk", m_xt4nRecTk );
00892 sc = m_xtupleEvt->addItem ("nTdsTk", m_xt4nTdsTk );
00893 sc = m_xtupleEvt->addItem ("t0", m_xt4t0);
00894 sc = m_xtupleEvt->addItem ("t0Stat", m_xt4t0Stat);
00895 sc = m_xtupleEvt->addItem ("t0Truth", m_xt4t0Truth);
00896 sc = m_xtupleEvt->addItem ("time", m_xt4time);
00897 sc = m_xtupleEvt->addItem ("nDigi", m_xt4nDigi, 0, 7000);//# of hit/mc track
00898 sc = m_xtupleEvt->addIndexedItem ("layer", m_xt4nDigi, m_xt4Layer);
00899 sc = m_xtupleEvt->addIndexedItem ("rt", m_xt4nDigi, m_xt4Time);
00900 sc = m_xtupleEvt->addIndexedItem ("rc", m_xt4nDigi, m_xt4Charge);
00901 //sc = m_xtupleEvt->addIndexedItem ("rawHit", m_xt4nDigi, m_xt4rawHit);
00902 //sc = m_xtupleEvt->addIndexedItem ("recHit", m_xt4nDigi, m_xt4recHit);
00903 } else { // did not manage to book the N tuple....
00904 log << MSG::ERROR << "Cannot book N-tuple: FILE121/evt" << endmsg;
00905 }
00906 }// end of book nt4
00907
00908 //book tuple of segment
00909 NTuplePtr ntSeg(ntupleSvc(), "FILE121/seg");
00910 if ( ntSeg ) { m_xtupleSeg = ntSeg;}
00911 else {
00912 m_xtupleSeg = ntupleSvc()->book ("FILE121/seg", CLID_ColumnWiseTuple, "MdcxTrackFinder segment data");
00913 if ( m_xtupleSeg ) {
00914 sc = m_xtupleSeg->addItem ("sl", m_xtsSl);
00915 sc = m_xtupleSeg->addItem ("d0", m_xtsD0);
00916 sc = m_xtupleSeg->addItem ("omega", m_xtsOmega);
00917 sc = m_xtupleSeg->addItem ("phi0", m_xtsPhi0);
00918 sc = m_xtupleSeg->addItem ("d0sl", m_xtsD0_sl_approx);
00919 sc = m_xtupleSeg->addItem ("phi0sl", m_xtsPhi0_sl_approx);
00920 sc = m_xtupleSeg->addItem ("xbbrrf", m_xtsXline_bbrrf);
00921 sc = m_xtupleSeg->addItem ("ybbrrf", m_xtsYline_bbrrf);
00922 sc = m_xtupleSeg->addItem ("xslope", m_xtsXline_slope);
00923 sc = m_xtupleSeg->addItem ("yslope", m_xtsYline_slope);
00924 sc = m_xtupleSeg->addItem ("pat", m_xtsPat);
00925 sc = m_xtupleSeg->addItem ("chisq", m_xtsChisq);
00926 sc = m_xtupleSeg->addItem ("nDigi", m_xtsNDigi, 0, 20);
00927 sc = m_xtupleSeg->addIndexedItem ("layer", m_xtsNDigi, m_xtsLayer);
00928 sc = m_xtupleSeg->addIndexedItem ("wire", m_xtsNDigi, m_xtsWire);
00929 sc = m_xtupleSeg->addIndexedItem ("inSeg", m_xtsNDigi, m_xtsInSeg);
00930 }
00931 }
00932 NTuplePtr nt5(ntupleSvc(), "FILE121/trkl");
00933 if ( nt5 ) { m_xtupleTrkl= nt5;}
00934 else {
00935 m_xtupleTrkl= ntupleSvc()->book ("FILE121/trkl", CLID_RowWiseTuple, "MdcxReco track info");
00936 if ( m_xtupleTrkl) {
00937 sc = m_xtupleTrkl->addItem ("layer", m_xt5Layer);
00938 sc = m_xtupleTrkl->addItem ("wire", m_xt5Wire);
00939 }
00940 }
00941
00942 NTuplePtr ntCsmcSew(ntupleSvc(), "FILE121/csmc");
00943 if ( ntCsmcSew ) { m_xtupleCsmcSew = ntCsmcSew;}
00944 else {
00945 m_xtupleCsmcSew = ntupleSvc()->book ("FILE121/csmc", CLID_ColumnWiseTuple, "MdcxReco reconsturction results");
00946 if ( m_xtupleCsmcSew ) {
00947 sc = m_xtupleCsmcSew->addItem ("dD0", m_csmcD0);
00948 sc = m_xtupleCsmcSew->addItem ("dPhi0", m_csmcPhi0);
00949 sc = m_xtupleCsmcSew->addItem ("dZ0", m_csmcZ0);
00950 sc = m_xtupleCsmcSew->addItem ("dOmega", m_csmcOmega);
00951 sc = m_xtupleCsmcSew->addItem ("dPt", m_csmcPt);
00952 sc = m_xtupleCsmcSew->addItem ("dTanl", m_csmcTanl);
00953 }
00954 }
00955 //--------------end of book ntuple------------------
00956 }
|
|
|
|
|
|
01253 {
01254 double tdr[43];
01255 double tdr_wire[43];
01256 for(int i=0; i<43; i++){tdr[i]=9999.;}
01257
01258 // make flag
01259 TrkHotList::hot_iterator hotIter = trkHitList->hotList().begin();
01260 while (hotIter!=trkHitList->hotList().end()) {
01261 MdcHitOnTrack* hot = const_cast<MdcHitOnTrack*> (&(*hotIter->mdcHitOnTrack()));
01262
01263 //driftTime(tof,z)
01264 double dt = hot->mdcHit()->driftTime(0.,0.);
01265 int layer = hot->mdcHit()->layernumber();
01266 int wire = hot->mdcHit()->wirenumber();
01267
01268 //std::cout<<__FILE__<<" "<<dt<< std::endl;
01269 if (dt < tdr[layer]) {
01270 tdr[layer] = dt;
01271 tdr_wire[layer] = wire;
01272 }
01273 hotIter++;
01274 }
01275
01276 //std::cout<<" tdr wire ";
01277 //for(int i=0;i<43;i++){
01278 //std::cout<<i<<","<<tdr_wire[i]<<" tdr="<<tdr[i]<<"\n";
01279 //}
01280 // inactive multi hit
01281 hotIter = trkHitList->hotList().begin();
01282 while (hotIter!=trkHitList->hotList().end()) {
01283 int layer = hotIter->mdcHitOnTrack()->mdcHit()->layernumber();
01284 int wire = hotIter->mdcHitOnTrack()->mdcHit()->wirenumber();
01285 //double dt = hotIter->mdcHitOnTrack()->mdcHit()->driftTime(0.,0.)-m_bunchT0;
01286
01287 if ((tdr[layer] <9998.) && (tdr_wire[layer]!=wire)){
01288 MdcHitOnTrack* hot = const_cast<MdcHitOnTrack*> (&(*hotIter->mdcHitOnTrack()));
01289 hot->setActivity(false);
01290 trkHitList->removeHit( hotIter->mdcHitOnTrack()->mdcHit() );
01291 //std::cout<<__FILE__<<" inactive "<< layer<<" "<<wire<<" dt "<<dt << std::endl;
01292 }else{
01293 hotIter++;
01294 }
01295 }
01296 }
|
|
|
|
|
|
01105 {
01106 cout << name()<<" found " <<segList.length() << " segs :" << endl;
01107 for (int i =0; i< segList.length(); i++){
01108 std::cout<<i<<" ";
01109 segList[i]->printSeg();
01110 }
01111 }
|
|
|
|
|
|
01346 {
01347 std::cout<<__FILE__<<" "<<__LINE__<<" All hits in TDS, nhit="<<hitList->size()<< std::endl;
01348 RecMdcHitCol::iterator it = hitList->begin();
01349 for(;it!= hitList->end(); it++){
01350 RecMdcHit* h = (*it);
01351 Identifier id(h->getMdcId());
01352 int layer = MdcID::layer(id);
01353 int wire = MdcID::wire(id);
01354 cout<<"("<< layer <<","<<wire<<") lr:"<<h->getFlagLR()<<" stat:"<<h->getStat()<<" tk:"<<h->getTrkId() <<" doca:"<<setw(10)<<h->getDoca()<<std::endl;
01355 }//end of hit list
01356 }
|
|
|
|
|
|
01297 {
01298 std::cout<< "tksize = "<<trackList->size() << std::endl;//yzhang debug
01299 RecMdcTrackCol::iterator it = trackList->begin();
01300 for (;it!= trackList->end();it++){
01301 RecMdcTrack *tk = *it;
01302 std::cout<< "//====RecMdcTrack "<<tk->trackId()<<"====:" << std::endl;
01303 cout <<" d0 "<<tk->helix(0)
01304 <<" phi0 "<<tk->helix(1)
01305 <<" cpa "<<tk->helix(2)
01306 <<" z0 "<<tk->helix(3)
01307 <<" tanl "<<tk->helix(4)
01308 <<endl;
01309 std::cout<<" q "<<tk->charge()
01310 <<" theta "<<tk->theta()
01311 <<" phi "<<tk->phi()
01312 <<" x0 "<<tk->x()
01313 <<" y0 "<<tk->y()
01314 <<" z0 "<<tk->z()
01315 <<" r0 "<<tk->r()
01316 <<endl;
01317 std::cout <<" p "<<tk->p()
01318 <<" pt "<<tk->pxy()
01319 <<" px "<<tk->px()
01320 <<" py "<<tk->py()
01321 <<" pz "<<tk->pz()
01322 <<endl;
01323 std::cout<<" tkStat "<<tk->stat()
01324 <<" chi2 "<<tk->chi2()
01325 <<" ndof "<<tk->ndof()
01326 <<" nhit "<<tk->getNhits()
01327 <<" nst "<<tk->nster()
01328 <<endl;
01329
01330 int nhits = tk->getVecHits().size();
01331 std::cout<<nhits <<" Hits: " << std::endl;
01332 for(int ii=0; ii <nhits ; ii++){
01333 Identifier id(tk->getVecHits()[ii]->getMdcId());
01334 int layer = MdcID::layer(id);
01335 int wire = MdcID::wire(id);
01336 cout<<"("<< layer <<","<<wire<<","<<tk->getVecHits()[ii]->getStat()
01337 <<",lr:"<<tk->getVecHits()[ii]->getFlagLR()<<") ";
01338 }//end of hit list
01339 std::cout << " "<< std::endl;
01340 }//end of tk list
01341 std::cout << " "<< std::endl;
01342
01343
01344 }//end of dumpTdsTrack
|
|
|
|
|
|
00726 {
00727 RecMdcTrackCol::iterator i_tk = trackList->begin();
00728 for (;i_tk != trackList->end(); i_tk++) {
00729 printTrack(*(i_tk));
00730 }
00731 }// dumpTrack
|
|
|
|
|
|
00242 {
00243
00244
00245 b_saveEvent=false;
00246 setFilterPassed(b_saveEvent);
00247 #ifdef MDCXTIMEDEBUG
00248 m_timer[0]->start();
00249 #endif
00250 MsgStream log(msgSvc(), name());
00251 log << MSG::INFO << "in execute()" << endreq;
00252 StatusCode sc;
00253
00254 nTk = 0; t_nTdsTk = 0; //yzhang for fill
00255 //------------------------------------
00256 // Get event No.
00257 //------------------------------------
00258 SmartDataPtr<Event::EventHeader> evtHead(eventSvc(),"/Event/EventHeader");
00259 if (!evtHead) {
00260 log << MSG::FATAL<< "Could not retrieve event header" << endreq;
00261 return StatusCode::FAILURE;
00262 }
00263 m_eventNo = evtHead->eventNumber();
00264 //std::cout << "x evt: " << m_eventNo<< std::endl;
00265 long t_evtNo = m_eventNo;
00266 //if (t_evtNo % 1000 == 0) std::cout << "x evt: " << t_evtNo << std::endl;
00267 IDataManagerSvc *dataManSvc;
00268 DataObject *aTrackCol;
00269 DataObject *aHitCol;
00270 if (!m_salvageTrk) {
00271 dataManSvc = dynamic_cast<IDataManagerSvc*> (eventSvc());
00272 eventSvc()->findObject("/Event/Recon/RecMdcTrackCol",aTrackCol);
00273 if(aTrackCol != NULL) {
00274 dataManSvc->clearSubTree("/Event/Recon/RecMdcTrackCol");
00275 eventSvc()->unregisterObject("/Event/Recon/RecMdcTrackCol");
00276 }
00277 eventSvc()->findObject("/Event/Recon/RecMdcHitCol",aHitCol);
00278 if(aHitCol != NULL) {
00279 dataManSvc->clearSubTree("/Event/Recon/RecMdcHitCol");
00280 eventSvc()->unregisterObject("/Event/Recon/RecMdcHitCol");
00281 }
00282 }
00283
00284 //------------------------------------
00285 // Initialize track collection in TDS
00286 //------------------------------------
00287 DataObject *aReconEvent;
00288 eventSvc()->findObject("/Event/Recon",aReconEvent);
00289 if (aReconEvent==NULL) {
00290 aReconEvent = new ReconEvent();
00291 sc = eventSvc()->registerObject("/Event/Recon",aReconEvent);
00292 if(sc != StatusCode::SUCCESS) {
00293 log << MSG::FATAL << "Could not register ReconEvent" <<endreq;
00294 return StatusCode::FAILURE;
00295 }
00296 }
00297 RecMdcTrackCol* trackList;
00298 eventSvc()->findObject("/Event/Recon/RecMdcTrackCol",aTrackCol);
00299 if (aTrackCol) {
00300 trackList = dynamic_cast<RecMdcTrackCol*> (aTrackCol);
00301 }else{
00302 trackList = new RecMdcTrackCol;
00303 sc = eventSvc()->registerObject(EventModel::Recon::RecMdcTrackCol, trackList);
00304 if(!sc.isSuccess()) {
00305 log << MSG::FATAL << " Could not register RecMdcTrack collection" <<endreq;
00306 return StatusCode::FAILURE;
00307 }
00308 }
00309 RecMdcHitCol* hitList;
00310 eventSvc()->findObject("/Event/Recon/RecMdcHitCol",aHitCol);
00311 if (aHitCol) {
00312 hitList = dynamic_cast<RecMdcHitCol*> (aHitCol);
00313 }else{
00314 hitList = new RecMdcHitCol;
00315 sc = eventSvc()->registerObject(EventModel::Recon::RecMdcHitCol, hitList);
00316 if(!sc.isSuccess()) {
00317 log << MSG::FATAL << " Could not register RecMdcHit collection" <<endreq;
00318 return StatusCode::FAILURE;
00319 }
00320 }
00321
00322 //------------------------------------
00323 // Initialize hit collection in TDS
00324 //------------------------------------
00325 DataObject *pnode = 0;
00326 sc = eventSvc()->retrieveObject("/Event/Hit",pnode);
00327 if(!sc.isSuccess()) {
00328 pnode = new DataObject;
00329 sc = eventSvc()->registerObject("/Event/Hit",pnode);
00330 if(!sc.isSuccess()) {
00331 log << MSG::FATAL << " Could not register /Event/Hit branch " <<endreq;
00332 return StatusCode::FAILURE;
00333 }
00334 }
00335 SmartDataPtr<MdcHitCol> m_hitCol(eventSvc(), "/Event/Hit/MdcHitCol");
00336 if (!m_hitCol){
00337 m_hitCol= new MdcHitCol;
00338 sc = eventSvc()->registerObject("/Event/Hit/MdcHitCol",m_hitCol);
00339 if(!sc.isSuccess()) {
00340 log << MSG::FATAL << " Could not register hit collection" <<endreq;
00341 return StatusCode::FAILURE;
00342 }
00343 }
00344
00345
00346 //------------------------------------
00347 // Get bunch time t0 (ns) and timing
00348 //------------------------------------
00349 m_bunchT0 = -999.;
00350 SmartDataPtr<RecEsTimeCol> aevtimeCol(eventSvc(),"/Event/Recon/RecEsTimeCol");
00351 if (!aevtimeCol || aevtimeCol->size()==0) {
00352 log << MSG::WARNING<< "evt "<<m_eventNo<<" Could not find RecEsTimeCol"<< endreq;
00353 return StatusCode::SUCCESS;
00354 }
00355
00356 RecEsTimeCol::iterator iter_evt = aevtimeCol->begin();
00357 for(; iter_evt!=aevtimeCol->end(); iter_evt++){
00358 m_bunchT0 = (*iter_evt)->getTest();
00359 m_t0Stat = (*iter_evt)->getStat();
00360 if ((m_t0Stat==0) || (m_bunchT0 < 0.) || (m_bunchT0 > 9999.0) ){
00361 log << MSG::WARNING << "Skip evt:"<<m_eventNo<< " by t0 = "<<m_bunchT0 << endreq;
00362 //return StatusCode::SUCCESS;
00363 }
00364 }
00365 int trigtiming=-10;
00366 SmartDataPtr<TrigData> trigData(eventSvc(),"/Event/Trig/TrigData");
00367 if(trigData){
00368 log << MSG::INFO <<"Trigger conditions 0--43:"<<endreq;
00369 for(int i = 0; i < 48; i++) {
00370 log << MSG::INFO << trigData->getTrigCondName(i)<<" ---- "<<trigData->getTrigCondition(i)<< endreq;
00371 }
00372 for(int i = 0; i < 16; i++) log << MSG::INFO << "Trigger channel "<< i << ": " << trigData->getTrigChannel(i) << endreq;
00373 m_timing=trigData->getTimingType();
00374 //cout<<"-----------------trigger timing type-----------------------: "<<trigtiming<<endl;
00375 log << MSG::INFO <<"Tigger Timing type: "<< trigtiming << endreq;
00376 }
00377
00378 //------------------------------------
00379 // Initialize MdcxHits
00380 //------------------------------------
00381 m_mdcxHits.reset();
00382 uint32_t getDigiFlag = 0;
00383 getDigiFlag += m_maxMdcDigi;
00384 if(m_dropHot||m_salvageTrk) getDigiFlag |= MdcRawDataProvider::b_dropHot;
00385 if(m_keepBadTdc) getDigiFlag |= MdcRawDataProvider::b_keepBadTdc;
00386 if(m_keepUnmatch) getDigiFlag |= MdcRawDataProvider::b_keepUnmatch;
00387 mdcDigiVec = m_rawDataProviderSvc->getMdcDigiVec(getDigiFlag);
00388 t_nDigi = mdcDigiVec.size();
00389
00390 if (0 == t_nDigi ){
00391 log << MSG::WARNING << " No hits in MdcDigiVec" << endreq;
00392 return StatusCode::SUCCESS;
00393 }
00394
00395
00396 // fill Mc truth
00397 //if(m_hist) fillMcTruth();
00398
00399 //skip event by hit numbe
00400 if (t_nDigi < m_minMdcDigi){
00401 log << MSG::WARNING << " Skip this event for MdcDigiVec.size() < "<<m_minMdcDigi << endreq;
00402 return StatusCode::SUCCESS;
00403 }
00404 m_mdcxHits.create(mdcDigiVec, m_bunchT0, m_cresol);
00405 const HepAList<MdcxHit>& dchitlist = m_mdcxHits.GetMdcxHitList();
00406
00407 //--------------------------------------------
00408 // Make segments (MdcxSeg's) out of MdcxHit's
00409 //--------------------------------------------
00410 MdcxFindSegs dcsegs(dchitlist,m_debug);
00411 const HepAList<MdcxSeg>& seglist = dcsegs.GetMdcxSeglist();
00412 if(m_debug > 1 ){ dumpMdcxSegs(seglist);}
00413 //if(m_hist){ fillMdcxSegs(seglist);}
00414
00415 //--------------------------------------------
00416 // Make tracks (MdcxFittedHel's) out of MdcxSeg's
00417 //--------------------------------------------
00418 MdcxFindTracks dctrks(seglist,m_debug);
00419 HepAList<MdcxFittedHel>& firsttrkl = (HepAList<MdcxFittedHel>&)dctrks.GetMdcxTrklist();
00420
00421 //if(m_hist){ fillTrkl(firsttrkl);}
00422
00423 MdcxMergeDups dcmergeem(firsttrkl);
00424 HepAList<MdcxFittedHel>& trkl = (HepAList<MdcxFittedHel>&)dcmergeem.GetMergedTrklist();
00425
00426 if (m_debug > 1 ) cout << "MdcxTrackFinder: after MergeDups, have "
00427 << trkl.length() << " track(s). nhits=" << dchitlist.length() << endl;
00428
00429 //---------------------------------------------------------
00430 // Put my tracks into official fitter and store to TDS
00431 //----------------------------------------------------------
00432
00433
00434 sc = FitMdcxTrack(trkl, dchitlist, m_hitCol, trackList, hitList);
00435 if (!sc.isSuccess()) {return StatusCode::SUCCESS;}
00436 t_nTdsTk = trackList->size();
00437
00438 #ifdef MDCXTIMEDEBUG
00439 m_timer[0]->stop();
00440 m_xt4time = m_timer[0]->elapsed();
00441 m_xtupleEvt->write();
00442 #endif
00443
00444 if(m_hist) fillEvent();
00445 if (m_debug > 0) {
00446 DataObject* pNode;
00447 eventSvc()->retrieveObject("/Event/Recon/RecMdcTrackCol",pNode);
00448 RecMdcTrackCol *tmpTrackCol = dynamic_cast<RecMdcTrackCol*> (pNode);
00449 eventSvc()->retrieveObject("/Event/Recon/RecMdcHitCol",pNode);
00450 int nTdsTk = 0;
00451 if(tmpTrackCol) nTdsTk = tmpTrackCol->size();
00452
00453 //if (t_evtNo % 1000 == 0) {
00454 std::cout<< "MdcxTrackFinder: evtNo "<< m_eventNo << " t0="<<m_bunchT0
00455 <<" Found " <<trkl.length()
00456 <<" keep "<< t_nTdsTk
00457 <<" finialy keep "<< nTdsTk;
00458
00459 int ndelete =0; trkl.length() - trackList->size();
00460 if( ndelete>0 ) std::cout <<" delete "<< ndelete;
00461 std::cout <<" track(s)" <<endl;
00462 //}
00463
00464 if(m_debug>1)dumpTdsTrack(tmpTrackCol);
00465 //dumpTdsHits(tmpHitCol);
00466
00467 }
00468 if((trackList->size()!=4) ) b_saveEvent = true;
00469 setFilterPassed(b_saveEvent);
00470 return StatusCode::SUCCESS;
00471 }
|
|
|
|
|
|
00959 {
00960 //-----------get raw digi-----------------------
00961 if (t_nDigi<=0) return;
00962 if (m_xtupleEvt == NULL ) return;
00963 m_xt4EvtNo = m_eventNo;
00964 m_xt4t0 = m_bunchT0;
00965 m_xt4t0Stat = m_t0Stat;
00966 m_xt4t0Truth = m_t0Truth;
00967 m_xt4nRecTk = nTk;
00968 m_xt4nTdsTk = t_nTdsTk;
00969 m_xt4nDigi = t_nDigi;
00970 int iDigi=0;
00971 MdcDigiCol::iterator iter = mdcDigiVec.begin();
00972 for (;iDigi<t_nDigi; iter++ ) {
00973 int l = MdcID::layer((*iter)->identify());
00974 m_xt4Layer[iDigi] = l;
00975 //int w = MdcID::wire((*iter)->identify());
00976 m_xt4Time[iDigi] = RawDataUtil::MdcTime((*iter)->getTimeChannel());
00977 m_xt4Charge[iDigi] = RawDataUtil::MdcCharge((*iter)->getChargeChannel());
00978 //m_xt4rawHit[l]++;
00979 iDigi++;
00980 }//end for iter
00981 m_xtupleEvt->write();
00982 }
|
|
|
|
|
|
01196 {
01197 MsgStream log(msgSvc(), name());
01198 StatusCode sc;
01199 //Initialize
01200 for (int ii=0;ii<43;ii++){
01201 for (int jj=0;jj<288;jj++){
01202 haveDigi[ii][jj] = -2;
01203 }
01204 }
01205 int nDigi = mdcDigiVec.size();
01206 for (int iDigi =0 ;iDigi<nDigi; iDigi++ ) {
01207 int l = MdcID::layer((mdcDigiVec[iDigi])->identify());
01208 int w = MdcID::wire((mdcDigiVec[iDigi])->identify());
01209 //haveDigi[l][w]=(mdcDigiVec[iDigi])->getTrackIndex();
01210 haveDigi[l][w]=1;
01211 }
01212
01213 if( m_mcHist ){
01214 //------------------get event start time truth-----------
01215 m_t0Truth = -10.;
01216 SmartDataPtr<Event::McParticleCol> mcParticleCol(eventSvc(),"/Event/MC/McParticleCol");
01217 if(!mcParticleCol){
01218 log << MSG::INFO << "Could not retrieve McParticelCol" << endreq;
01219 }else {
01220 Event::McParticleCol::iterator iter_mc = mcParticleCol->begin();
01221 for (;iter_mc != mcParticleCol->end(); iter_mc++){
01222 if ((*iter_mc)->primaryParticle()){
01223 m_t0Truth = (*iter_mc)->initialPosition().t();
01224 //px = (*iter_mc)->initialFourMomentum().x()/1000.;//GeV
01225 //py = (*iter_mc)->initialFourMomentum().y()/1000.;//GeV
01226 //pz = (*iter_mc)->initialFourMomentum().z()/1000.;//GeV
01227 }
01228 }
01229 }
01230 }
01231 //------------------Retrieve MC truth MdcMcHit------------
01232 /*
01233 SmartDataPtr<Event::MdcMcHitCol> mcMdcMcHitCol(eventSvc(),"/Event/MC/MdcMcHitCol");
01234 if (!mcMdcMcHitCol) {
01235 log << MSG::WARNING << "Could not find MdcMcHit" << endreq;
01236 }else{
01237 Event::MdcMcHitCol::iterator iter_mchit = mcMdcMcHitCol->begin();
01238 for (;iter_mchit != mcMdcMcHitCol->end(); iter_mchit++ ) {
01239 const Identifier id= (*iter_mchit)->identify();
01240 int layer = MdcID::layer(id);
01241 int wire = MdcID::wire(id);
01242 mcDrift[layer][wire] = (*iter_mchit)->getDriftDistance(); //drift in MC.
01243 mcLR[layer][wire] = (*iter_mchit)->getPositionFlag();
01244 mcX[layer][wire] = (*iter_mchit)->getPositionX();
01245 mcY[layer][wire] = (*iter_mchit)->getPositionY();
01246 mcZ[layer][wire] = (*iter_mchit)->getPositionZ();
01247 if (mcLR[layer][wire] == 0) mcLR[layer][wire] = -1;
01248 }
01249 }
01250 */
01251 }
|
|
|
|
|
|
01113 {
01114
01115 int cellMax[43] ={
01116 40,44,48,56, 64,72,80,80, 76,76,88,88,
01117 100,100,112,112, 128,128,140,140, 160,160,160,160,
01118 176,176,176,176, 208,208,208,208, 240,240,240,240,
01119 256,256,256,256, 288,288,288 };
01120 // Fill hits of every layer after segment finding
01121 int hitInSegList[43][288];
01122 for (int ii=0;ii<43;ii++){
01123 for (int jj=0;jj<cellMax[ii];jj++){ hitInSegList[ii][jj] = 0; }
01124 }
01125 for (int i = 0; i < segList.length(); i++) {
01126 MdcxSeg* aSeg = segList[i];
01127 for (int ihit=0 ; ihit< aSeg->Nhits() ; ihit++){
01128 const MdcxHit* hit = aSeg->XHitList()[ihit];
01129 int layer = hit->Layer();
01130 int wire = hit->WireNo();
01131 hitInSegList[layer][wire]++;
01132 }
01133 }
01134 for (int i = 0; i < segList.length(); i++) {
01135 MdcxSeg* aSeg = segList[i];
01136 if(aSeg==NULL)continue;
01137 m_xtsSl = aSeg->SuperLayer();
01138 m_xtsD0 = aSeg->D0();
01139 m_xtsOmega = aSeg->Omega();
01140 m_xtsPhi0 = aSeg->Phi0();
01141 m_xtsD0_sl_approx = aSeg->D0_sl_approx();
01142 m_xtsPhi0_sl_approx = aSeg->Phi0_sl_approx();
01143 m_xtsXline_bbrrf = aSeg->Xline_bbrrf();
01144 m_xtsYline_bbrrf = aSeg->Yline_bbrrf();
01145 m_xtsXline_slope = aSeg->Xline_slope();
01146 m_xtsYline_slope = aSeg->Yline_slope();
01147 int patIndex = -1;
01148 for (int ii = 0;ii<14;ii++){
01149 if (aSeg->Pat() == (int) MdcxSegPatterns::patt4[ii]){
01150 patIndex = ii;
01151 break;
01152 }
01153 }
01154 for (int ii = 0;ii<20;ii++){
01155 if (aSeg->Pat() == (int) MdcxSegPatterns::patt3[ii]){
01156 patIndex = ii +15;
01157 break;
01158 }
01159 }
01160 m_xtsPat = patIndex;
01161 m_xtsChisq = aSeg->Chisq();
01162 m_xtsNDigi = aSeg->Nhits();
01163 for (int ihit=0 ; ihit< aSeg->Nhits() ; ihit++){
01164 const MdcxHit* hit = aSeg->XHitList()[ihit];
01165 int layer = hit->Layer();
01166 int wire = hit->WireNo();
01167 m_xtsLayer[ihit] = layer;
01168 m_xtsWire[ihit] = wire;
01169 m_xtsInSeg[ihit] = hitInSegList[layer][wire];
01170 //std::cout<< "hitInSegList "<<hitInSegList[layer][wire] << std::endl;//yzhang debug
01171 }
01172 m_xtupleSeg->write();
01173 }
01174
01175 }//end of fillMdcxSegs
|
|
|
|
|
|
00984 {
00985
00986 //-----------get MC truth data-------------------
00987 m_xevtNo = m_eventNo;
00988 int recHitMap[43][288]={0};
00989 //int nTk = (*m_tracks).nTrack();
00990 //for (int itrack = 0; itrack < nTk; itrack++) {
00991 if (atrack==NULL) return;
00992
00993 const TrkFit* fitResult = atrack->fitResult();
00994 if (fitResult == 0) return;//check the fit worked
00995
00996 //for fill ntuples
00997 int nSt=0; //int nSeg=0;
00998 int seg[11] = {0}; int segme;
00999 //-----------hit list-------------
01000 HitRefVec hitRefVec;
01001 const TrkHitList* hitlist = atrack->hits();
01002
01003 TrkHitList::hot_iterator hot = hitlist->begin();
01004 int layerCount[43]={0};
01005 int i=0;
01006 for (;hot!= hitlist->end();hot++){
01007
01008 const MdcRecoHitOnTrack* recoHot;
01009 recoHot = dynamic_cast<const MdcRecoHitOnTrack*> (&(*hot));
01010 int layer = recoHot->mdcHit()->layernumber();
01011 int wire = recoHot->mdcHit()->wirenumber();
01012 m_xlayer[i] = layer;
01013 //m_xt4recHit[layer]++;//fill rec hit for hit eff
01014 m_xwire[i] = wire;
01015 layerCount[layer]++;
01016 recHitMap[layer][wire]++;
01017 m_xambig[i] = recoHot->wireAmbig();// wire ambig
01018 //fltLen
01019 double fltLen = recoHot->fltLen();
01020 m_xfltLen[i] = fltLen;
01021 double tof = recoHot->getParentRep()->arrivalTime(fltLen);
01022 //position
01023 HepPoint3D pos; Hep3Vector dir;
01024 recoHot->getParentRep()->traj().getInfo(fltLen,pos,dir);
01025 m_xx[i] = pos.x();
01026 m_xy[i] = pos.y();
01027 m_xz[i] = pos.z();
01028 m_xdriftT[i] = recoHot->mdcHit()->driftTime(tof,pos.z());
01029 m_xtof[i] = tof;
01030 m_xdriftD[i] = recoHot->drift();
01031 m_xsigma[i] = recoHot->hitRms();
01032 m_xtdc[i] = recoHot->rawTime();
01033 m_xadc[i] = recoHot->mdcHit()->charge();
01034 m_xdoca[i] = recoHot->dcaToWire();//sign w.r.t. dirft() FIXME
01035 m_xentra[i] = recoHot->entranceAngle();
01036 //m_xentraHit[i] = recoHot->entranceAngleHit();
01037 m_xact[i] = recoHot->isActive();
01038 //resid
01039 double res=999.,rese=999.;
01040 if (recoHot->resid(res,rese,false)){
01041 }else{}
01042 m_xresid[i] = res;
01043 //for n seg
01044 segme=0;
01045 if ( layer >0 ) {segme=(layer-1)/4;}
01046 seg[segme]++;
01047 if (recoHot->layer()->view()) { ++nSt; }
01048 i++;
01049 }// end fill hit
01050
01051 int nSlay=0;
01052 for(int i=0;i<11;i++){
01053 if (seg[i]>0) nSlay++;
01054 }
01055 m_xnSlay = nSlay;
01056
01057 //------------fill track result-------------
01058 //m_xtkId = itrack;
01059 //track parameters at closest approach to beamline
01060 double fltLenPoca = 0.0;
01061 TrkExchangePar helix = fitResult->helix(fltLenPoca);
01062 m_xphi0 = helix.phi0();
01063 m_xtanl = helix.tanDip();
01064 m_xz0 = helix.z0();
01065 m_xd0 = helix.d0();
01066 if(fabs(m_xz0)>20||fabs(m_xd0)>2) {
01067 //b_saveEvent = true;
01068 if(m_debug>1) std::cout<<"evt:"<<m_eventNo<<" BigVtx:"
01069 <<" d0 "<<helix.d0()<<" z0 "<<helix.z0()<<std::endl;
01070 }
01071 m_xpt = fitResult->pt();
01072 if (m_xpt > 0.00001){
01073 m_xcpa = fitResult->charge()/fitResult->pt();
01074 }
01075 //momenta and position
01076 CLHEP::Hep3Vector p1 = fitResult->momentum(fltLenPoca);
01077 double px,py,pz,pxy;
01078 pxy = fitResult->pt();
01079 px = p1.x();
01080 py = p1.y();
01081 pz = p1.z();
01082 m_xp = p1.mag();
01083 m_xpz = pz;
01084 HepPoint3D poca = fitResult->position(fltLenPoca);
01085 m_xpocax = poca.x();
01086 m_xpocay = poca.y();
01087 m_xpocaz = poca.z();
01088
01089 m_xq = fitResult->charge();
01090 m_xnAct = fitResult->nActive();
01091 m_xnHit = hitlist->nHit();
01092 m_xnDof = fitResult->nDof();
01093 m_xnSt = nSt;
01094 m_xchi2 = fitResult->chisq();
01095 //for (int l=0;l<43;l++) m_xlayerCount[l] = layerCount[l];
01096 m_xt0 = m_bunchT0;
01097 m_xtiming = m_timing;
01098 m_xt0Stat = m_t0Stat;
01099 m_xt0Truth= m_t0Truth;
01100 m_xtuple1->write();
01101 //}//end of loop rec trk list
01102
01103 }//end of dumpTrackList
|
|
|
|
|
|
01177 {
01178 int nDigi = 0;
01179 int iDigi = 0;
01180 for (int i =0; i< firsttrkl.length(); i++){
01181 nDigi+= firsttrkl[i]->Nhits();
01182 }
01183 for (int i=0;i<firsttrkl.length();i++){
01184 for (int ihit=0 ; ihit< firsttrkl[i]->Nhits() ; ihit++){
01185 const MdcxHit* hit = firsttrkl[i]->XHitList()[ihit];
01186 int layer = hit->Layer();
01187 int wire = hit->WireNo();
01188 m_xt5Layer = layer;
01189 m_xt5Wire = wire;
01190 m_xtupleTrkl->write();
01191 iDigi++;
01192 }
01193 }
01194 }//end of fillTrkl
|
|
|
|
|
|
00473 {
00474 MsgStream log(msgSvc(), name());
00475 log << MSG::INFO << "in finalize()" << endreq;
00476 //tot evtNo, trkNum
00477 return StatusCode::SUCCESS;
00478 }
|
|
||||||||||||||||||||||||
|
|
|
||||||||||||||||||||||||
|
m_bunchT0 in "ns" here, but "second" in TrkLineMaker&TrkHelixMaker 00516 {
00517 StatusCode sc;
00518 MsgStream log(msgSvc(), name());
00519
00520 //--Add Hit to MdcxFittedHel
00521 MdcxAddHits dcaddem(trkl, dchitlist, m_addHitCut);
00522 if (m_debug > 1) cout << "MdcxTrackFinder: after AddHits, have "
00523 << trkl.length() << " track(s). nhits=" << dchitlist.length() << endl;
00524
00525 TrkLineMaker linefactory;
00526 TrkHelixMaker helixfactory;
00527
00528
00529 int tkLen = trkl.length();//FIXME
00530 for (int kk=0; kk< tkLen; kk++){
00531 const HepAList<MdcxHit>& xhits = trkl[kk]->XHitList();
00532 if(m_debug>2){
00533 std::cout<<__FILE__<<" FitMdcxTrack "<<kk<< std::endl;
00534 for(int i=0; i<xhits.length(); i++){ xhits[i]->print(std::cout); }
00535 }
00536
00537 if(m_debug>2) std::cout<<__FILE__<<" before add hits nhits="<<xhits.length()<< std::endl;
00538 HepAList<MdcxHit> xass = dcaddem.GetAssociates(kk);
00539 if(m_debug>2) std::cout<<__FILE__<<" after,add "<<xass.length()<<" hits"<<std::endl;
00540 MdcxFittedHel mdcxHelix = *trkl[kk];
00541 double thechisq = mdcxHelix.Chisq();
00542 TrkExchangePar tt(-mdcxHelix.D0(),mdcxHelix.Phi0(),mdcxHelix.Omega(),-mdcxHelix.Z0(),-mdcxHelix.Tanl());
00543 TrkRecoTrk* aTrack;
00544 int nparm;
00545
00546 if (m_lineFit) {
00548 aTrack = linefactory.makeTrack(tt, thechisq, *m_context, m_bunchT0*1.e-9);
00549 nparm = 4;
00550 linefactory.setFlipAndDrop(*aTrack, true, true);
00551 } else {
00552 aTrack = helixfactory.makeTrack(tt, thechisq, *m_context, m_bunchT0*1.e-9);
00553 nparm = 5;
00554 helixfactory.setFlipAndDrop(*aTrack, true, true);
00555 }
00556
00557 TrkHitList* m_trkHitList = aTrack->hits();
00558 if (0 == m_trkHitList) {
00559 delete aTrack;
00560 aTrack = NULL;
00561 continue;
00562 }
00563
00564
00565 MdcxHitsToHots(mdcxHelix, xhits, m_trkHitList, m_hitCol);
00566 //std::cout<<"xhits---------------------"<<std::endl;//debug
00567 //m_trkHitList->hotList().printAll(cout);//debug
00568 MdcxHitsToHots(mdcxHelix, xass, m_trkHitList, m_hitCol);
00569 //std::cout<<"xass----------------------"<<std::endl;//debug
00570 //m_trkHitList->hotList().printAll(cout);//debug
00571 //std::cout<<__FILE__<<" "<<__LINE__<< std::endl;
00572 //std::cout<<"size "<<m_trkHitList->hotList().nHit()<<std::endl;
00573 //int beforDrop = m_trkHitList->hotList().nHit();
00574 if(m_dropMultiHotInLayer) dropMultiHotInLayer(m_trkHitList);//yzhang debug FIXME
00575 //int afterDrop = m_trkHitList->hotList().nHit();
00576 //std::cout<<"drop "<<beforDrop-afterDrop<<" keep:"<<afterDrop<<::endl;
00577 TrkErrCode err = m_trkHitList->fit();
00578
00579 const TrkFit* theFit = aTrack->fitResult();
00580 float rcs = 10000.0;
00581
00582
00583 if (theFit) {
00584 int ndof = theFit->nActive()-nparm;
00585 if (ndof > 0) rcs = theFit->chisq()/float(ndof);
00586 if (m_debug>1) {
00587 if (4 == nparm) cout << " TrkLineMaker";
00588 else cout << " TrkHelixMaker";
00589 cout << " success trkNo. " << kk << " status " << err.success() << " rcs " << rcs
00590 << " chi2 " << theFit->chisq() << " nactive " << theFit->nActive() << endl;
00591 }
00592 }
00593 if ( (1 == err.success()) && (rcs < 20.0) ) {
00594 if(m_debug>1) std::cout<<"aTrack->fitResult() success "<<std::endl;//yzhang debug
00595 if (4 == nparm) {
00596 linefactory.setFlipAndDrop(*aTrack, false, false);
00597 } else {
00598 helixfactory.setFlipAndDrop(*aTrack, false, false);
00599 }
00600 //-------------Stick the found tracks into the list in RecMdcTrackCol--------
00601 if (m_debug>1) { cout << "MdcxTrackFinder: accept a track " << endl; }
00602 // update history
00603 aTrack->status()->addHistory(err,name().c_str()); //yzhang FIXME
00604 store(aTrack, trackList, hitList);//aTrack have been deleted
00605 } else if ( (2 == err.success()) && (rcs < 150.0) ) {
00606 if(m_debug > 1) std::cout<<"info:( err.success = 2, refit now)"<<std::endl;//yzhang debug
00607 int nrefit = 0;
00608 while (nrefit++ < 5) {
00609 if (m_debug>1) std::cout << "refit time " << nrefit << std::endl;
00610 err = m_trkHitList->fit();
00611 if (err.success() == 1) break;
00612 }
00613 if (err.success() == 1) {
00614 if (4 == nparm) {
00615 linefactory.setFlipAndDrop(*aTrack, false, false);
00616 } else {
00617 helixfactory.setFlipAndDrop(*aTrack, false, false);
00618 }
00619 //-------------Stick the found tracks into the list in RecMdcTrackCol--------
00620 //if (m_debug>1) { cout << "MdcxTrackFinder: accept a track and store to TDS" << endl; }
00621 // update history
00622 aTrack->status()->addHistory(err,name().c_str()); //yzhang FIXME
00623 store(aTrack, trackList, hitList);//aTrack have been deleted
00624 }
00625 } else {
00626 if (m_debug >1) {
00627 std::cout<<"info:( aTrack->fitResult() faild )"<<std::endl;//yzhang debug
00628 err.print(cout);
00629 cout << endl;
00630 }
00631 delete aTrack;
00632 aTrack = NULL;
00633 //---------------------------------------
00634 // Fit no good; try a better input helix
00635 //---------------------------------------
00636 if(m_debug>1) std::cout<<"info:( aTrack->fitResult() Fit no good; try a better input helix)"<<std::endl;//yzhang debug
00637 mdcxHelix.Grow(*trkl[kk],xass);
00638 mdcxHelix.VaryRes();
00639 mdcxHelix.SetChiDofBail(1500);//yzhang add 2009-11-03
00640 int fail = mdcxHelix.ReFit();
00641 if(m_debug>1)std::cout<<__FILE__<<" refit fail:"<<fail<< std::endl;
00642 if (!mdcxHelix.Fail()) {
00643 const HepAList<MdcxHit>& bxhits = mdcxHelix.XHitList();
00644 thechisq = mdcxHelix.Chisq();
00645 TrkExchangePar tb(mdcxHelix.D0(),mdcxHelix.Phi0(),mdcxHelix.Omega(),mdcxHelix.Z0(),mdcxHelix.Tanl());
00646 TrkRecoTrk* bTrack;
00647 if (4 == nparm){
00648 bTrack = linefactory.makeTrack(tb,thechisq,*m_context,m_bunchT0*1.e-9);
00649 linefactory.setFlipAndDrop(*bTrack, false, false);
00650 }else{
00651 bTrack = helixfactory.makeTrack(tb,thechisq,*m_context,m_bunchT0*1.e-9);
00652 helixfactory.setFlipAndDrop(*bTrack, false, false);
00653 }
00654 TrkHitList* bhits = bTrack->hits();
00655 if (0 == bhits){delete bTrack; bTrack = NULL; continue;}
00656
00657 MdcxHitsToHots(mdcxHelix, bxhits, bhits, m_hitCol);
00658 TrkErrCode berr = bhits->fit();
00659 const TrkFit* bFit = bTrack->fitResult();
00660 rcs=10000.0;
00661 if (bFit) {
00662 int ndof = bFit->nActive() - nparm;
00663 if (ndof > 0) rcs = bFit->chisq()/float(ndof);
00664 if (m_debug >1) {
00665 if (4 == nparm) cout << " TrkLineMaker";
00666 else cout << " TrkHelixMaker";
00667 cout << " success trkNo. " << kk << " status " << berr.success() << " rcs " << rcs
00668 << " chi2 " << bFit->chisq() << " nactive "<< bFit->nActive() << endl;
00669 }
00670 }
00671 if ( ( 1 == berr.success() ) && ( rcs < 50.0 ) ) {
00672 // update history
00673 bTrack->status()->addHistory(berr,name().c_str());//yzhang FIXME
00674 if (m_debug>1) {
00675 cout << "MdcxTrackFinder: accept b track and store to TDS" << endl;
00676 bTrack->printAll(cout);
00677 }
00678 store(bTrack, trackList, hitList);//bTrack have been deleted
00679 } else {
00680 if (m_debug>1) {
00681 cout<< " fit failed "<<endl;
00682 berr.print(cout);
00683 cout << endl;
00684 }
00685 if (bTrack!=NULL) { delete bTrack; bTrack = NULL; }
00686 }
00687 }else{
00688 //cout<< " grow and refit failed "<<endl;
00689 }
00690 }
00691 }
00692 if(m_debug >1) dumpTdsTrack(trackList);
00693 return StatusCode::SUCCESS;
00694
00695 }// end of FitMdcxTrack
|
|
|
|
|
|
00165 {
00166 MsgStream log(msgSvc(), name());
00167 log << MSG::INFO << "in initialize()" << endreq;
00168
00169
00170 //m_flags.readPar(m_paramFile);
00171 #ifdef MDCXTIMEDEBUG
00172 StatusCode tsc = service( "BesTimerSvc", m_timersvc);
00173 if( tsc.isFailure() ) {
00174 log << MSG::WARNING << name() << ": Unable to locate BesTimer Service" << endreq;
00175 return StatusCode::FAILURE;
00176 }
00177 m_timer[0] = m_timersvc->addItem("Execution");
00178 m_timer[0]->propName("Execution");
00179 #endif
00180
00181 if (m_helixFitCut.size() == 43){
00182 for(int i=0; i<43; i++){
00183 MdcTrkReconCut_helix_fit[i] = m_helixFitCut[i];
00184 }
00185 }else{
00186 for(int i=0; i<43; i++){
00187 MdcTrkReconCut_helix_fit[i] = 5.;
00188 }
00189 }
00190 MdcxParameters::debug = m_debug;
00191 MdcxParameters::minTrkProb = m_minTrkProb;
00192 MdcxParameters::csmax4 = m_csmax4;
00193 MdcxParameters::csmax3 = m_csmax3;
00194 MdcxParameters::helixFitSigma = m_helixFitSigma;
00195 MdcxParameters::maxRcsInAddSeg= m_maxRcsInAddSeg;
00196 MdcxParameters::nSigAddHitTrk = m_nSigAddHitTrk;
00197 MdcxParameters::maxProca = m_maxProca;
00198 Pdt::readMCppTable(m_pdtFile);
00199 MdcxFittedHel::debug = m_debug;
00200
00201
00202 // Get MdcCalibFunSvc
00203 IMdcCalibFunSvc* imdcCalibSvc;
00204 StatusCode sc = service ("MdcCalibFunSvc", imdcCalibSvc);
00205 m_mdcCalibFunSvc = dynamic_cast<MdcCalibFunSvc*>(imdcCalibSvc);
00206 if ( sc.isFailure() ){
00207 log << MSG::FATAL << "Could not load MdcCalibFunSvc!" << endreq;
00208 return StatusCode::FAILURE;
00209 }
00210 MdcxHit::setMdcCalibFunSvc(m_mdcCalibFunSvc);
00211 MdcxHit::setCountPropTime(m_countPropTime);
00212
00213 // Get RawDataProviderSvc
00214 IRawDataProviderSvc* iRawDataProvider;
00215 sc = service ("RawDataProviderSvc", iRawDataProvider);
00216 if ( sc.isFailure() ){
00217 log << MSG::FATAL << "Could not load RawDataProviderSvc!" << endreq;
00218 return StatusCode::FAILURE;
00219 }
00220 m_rawDataProviderSvc = dynamic_cast<RawDataProviderSvc*>(iRawDataProvider);
00221
00222
00223 //Initailize magnetic filed
00224 sc = service ("MagneticFieldSvc",m_pIMF);
00225 if(sc!=StatusCode::SUCCESS) {
00226 log << MSG::ERROR << "Unable to open Magnetic field service"<<endreq;
00227 }
00228 m_bfield = new BField(m_pIMF);
00229 log << MSG::INFO << "field z = "<<m_bfield->bFieldNominal()<< endreq;
00230 m_context = new TrkContextEv(m_bfield);
00231
00232 if (m_hist) {bookNTuple();}
00233 if (m_dropHitsSigma.size()==43){
00234 for (int ii=0;ii<43;ii++) {
00235 MdcxParameters::dropHitsSigma[ii]=m_dropHitsSigma[ii];
00236 }
00237 }
00238
00239 return StatusCode::SUCCESS;
00240 }
|
|
||||||||||||||||
|
|
|
||||||||||||||||||||
|
|
|
||||||||||||||||
|
|
|
||||||||||||||||||||
|
00482 {
00483
00484 if ( 0 == mdcxHits.length() ) return;
00485
00486 int ihits = 0;
00487 while(mdcxHits[ihits]) {
00488 int ambig=mdcxHelix.Doca_Samb();
00489 const MdcHit* newhit = mdcxHits[ihits]->getMdcHit();
00490 if ( 0 == newhit ) {
00491 const MdcDigi* theDigi = mdcxHits[ihits]->getDigi();
00492 int layer = MdcID::layer(mdcxHits[ihits]->getDigi()->identify());
00493 int wire = MdcID::wire(mdcxHits[ihits]->getDigi()->identify());
00494 m_digiMap[layer][wire] = mdcxHits[ihits]->getDigi();
00495 MdcHit *thehit = new MdcHit(theDigi, m_gm);
00496 thehit->setCalibSvc(m_mdcCalibFunSvc);
00497 thehit->setCountPropTime(m_countPropTime);
00498 thehit->setCosmicFit(m_cosmicFit);
00499
00500 mdcHitCol->push_back(thehit);
00501 newhit = thehit;
00502 }
00503 MdcRecoHitOnTrack temp(*newhit, ambig, m_bunchT0);//m_bunchT0 nano second here
00504 MdcHitOnTrack* newhot = &temp;
00505 double fltLen = mdcxHelix.Doca_FLen();
00506 newhot->setFltLen(fltLen);
00507
00508 // Store MdcxHits to TrkHitList
00509 m_trkHitList->appendHot(newhot); //yzhang TEMP FIXME
00510 ihits++;
00511 }
00512 }
|
|
|
|
|
|
00733 {
00734 //yzhang debug
00735 std::cout<< " MdcTrack Id:"<<tk->trackId() <<" q:"<< tk->charge()<< std::endl;
00736 std::cout<< "dr Fi0 Cpa Dz Tanl Chi2 Ndf nSt FiTerm poca" << std::endl;
00737 std::cout<<"(" <<setw(5) << tk->helix(0)<<","<< setw(5) << tk->helix(1)<<"," << setw(5) << tk->helix(2) <<","
00738 << setw(5) << tk->helix(3) << ","<< setw(5) << tk->helix(4) <<")"
00739 << setw(5) << tk->chi2() << setw(4) << tk->ndof()
00740 << setw(4) << tk->getNhits() << setw(4) << tk->nster()
00741 << setw(5) << tk->getFiTerm() <<tk->poca()<<std::endl;
00742 std::cout<< " ErrMat "<<tk->err() << std::endl;
00743
00744 std::cout<< "hitId tkId (l,w) fltLen lr dt ddl ddre tdc chi2Add doca entr z tprop stat " << std::endl;
00745
00746 HitRefVec hl = tk->getVecHits();
00747 HitRefVec::iterator it = hl.begin();
00748 for (;it!=hl.end();++it){
00749 RecMdcHit* h = *it;
00750 int layer = MdcID::layer(h->getMdcId());
00751 double _vprop = (layer<8) ? Constants::vpropInner : Constants::vpropOuter;
00752 const MdcLayer* _layerPtr = m_gm->Layer(layer);
00753 double _zlen = _layerPtr->zLength();
00754 double z = h->getZhit();
00755 double tprop;
00756 if (0 == layer%2){
00757 tprop = (0.5*_zlen + z)/_vprop; //odd
00758 }else{
00759 tprop = (0.5*_zlen - z)/_vprop; //even
00760 }
00761 // build the sense wires
00762 //const MdcSWire* wire = m_gm->Wire(MdcID::layer(h->getMdcId()),MdcID::wire(h->getMdcId()));
00763 //double z = wire->zForw();
00764 //while(z<wire->zRear()){
00765 //HepPoint3D pos;
00766 //Hep3Vector dir;
00767 //if(!(wire->getTraj()==NULL)){
00768 //wire->getTraj()->getInfo(z,pos,dir);
00769 //}
00770 //std::cout<<"("<< wire->layer()->layNum()<<","<<wire->cell()<<" "<<wire->Id()<<")";
00771 //std::cout<<" z, sag:"<<z
00772 //<<", "<<wire->getTraj()->deltaY(z-wire->zForw())<<std::endl;
00774 //z+=1.;
00775 //}
00776 std::cout<< setw(4) << h->getId() << setw(4) << h->getTrkId() <<
00777 setw(4) << MdcID::layer(h->getMdcId()) <<setw(5) << MdcID::wire(h->getMdcId()) <<
00778 setw(10) << h->getFltLen() <<
00779 setw(3) << h->getFlagLR() <<setw(10) << h->getDriftT() <<
00780 setw(12) << h->getDriftDistLeft() <<setw(8) << h->getErrDriftDistRight() <<
00781 setw(8) << h->getTdc() <<setw(8) << h->getChisqAdd() <<
00782 setw(10) << h->getDoca() <<setw(10) << h->getEntra() <<
00783 setw(10) << h->getZhit() << setw(10) << tprop<<
00784 setw(5)<< h->getStat() << std::endl;
00785 }
00786 }//print track
|
|
||||||||||||||||
|
|
|
||||||||||||||||
|
00698 {
00699 assert (aTrack != NULL);
00700 nTk++;
00701 int trackId = trackList->size();
00702 TrkExchangePar helix = aTrack->fitResult()->helix(0.);
00703 if(m_dropTrkPt>0. && (aTrack->fitResult()->pt()<m_dropTrkPt)) {
00704 //std::cout<< __FILE__ <<" delete track by pt "
00705 //<<aTrack->fitResult()->pt()<<"<ptCut "<<m_dropTrkPt << std::endl;
00706 return;
00707 }
00708
00709 if( ( (fabs(helix.d0())>m_d0Cut) ||( fabs(helix.z0())>m_z0Cut) ) ){
00710 //std::cout<< __FILE__ <<" delete track by d0 "<<helix.d0()<<">d0Cut "<<m_d0Cut
00711 //<<" or z0 "<<helix.z0()<<" >z0Cut "<<m_z0Cut << std::endl;
00712 return;
00713 }
00714
00715 if(m_hist) fillTrack(aTrack);
00716 MdcTrack mdcTrack(aTrack);//aTrack have been deleted in ~MdcTrack()
00717 //tkStat: 0,PatRec 1,MdcxReco 2,Tsf 3,CurlFinder -1,Combined cosmic
00718 int tkStat = 1;
00719 int nHitbefore = hitList->size();
00720
00721 mdcTrack.storeTrack(trackId, trackList, hitList, tkStat);
00722 int nHitAfter = hitList->size();
00723 if (nHitAfter - nHitbefore <10 ) setFilterPassed(true);
00724 }
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
1.3.9.1