346 NTuplePtr nt1(
ntupleSvc(),
"FILE104/n101");
348 if ( nt1 ) m_nt1 = nt1;
350 m_nt1=
ntupleSvc()->book(
"FILE104/n101",CLID_ColumnWiseTuple,
"KalFit");
353 status = m_nt1->addItem(
"trackid",m_trackid);
354 status = m_nt1->addItem(
"stat",5,2,m_stat);
355 status = m_nt1->addItem(
"ndf",5,2,m_ndf);
356 status = m_nt1->addItem(
"chisq",5,2,m_chisq);
357 status = m_nt1->addItem(
"length",5,m_length);
358 status = m_nt1->addItem(
"tof",5,m_tof);
359 status = m_nt1->addItem(
"nhits",5,m_nhits);
360 status = m_nt1->addItem(
"zhelix",5,m_zhelix);
361 status = m_nt1->addItem(
"zhelixe",5,m_zhelixe);
362 status = m_nt1->addItem(
"zhelixmu",5,m_zhelixmu);
363 status = m_nt1->addItem(
"zhelixk",5,m_zhelixk);
364 status = m_nt1->addItem(
"zhelixp",5,m_zhelixp);
365 status = m_nt1->addItem(
"zptot",m_zptot);
366 status = m_nt1->addItem(
"zptote",m_zptote);
367 status = m_nt1->addItem(
"zptotmu",m_zptotmu);
368 status = m_nt1->addItem(
"zptotk",m_zptotk);
369 status = m_nt1->addItem(
"zptotp",m_zptotp);
371 status = m_nt1->addItem(
"zpt",m_zpt);
372 status = m_nt1->addItem(
"zpte",m_zpte);
373 status = m_nt1->addItem(
"zptmu",m_zptmu);
374 status = m_nt1->addItem(
"zptk",m_zptk);
375 status = m_nt1->addItem(
"zptp",m_zptp);
377 status = m_nt1->addItem(
"fptot",m_fptot);
378 status = m_nt1->addItem(
"fptote",m_fptote);
379 status = m_nt1->addItem(
"fptotmu",m_fptotmu);
380 status = m_nt1->addItem(
"fptotk",m_fptotk);
381 status = m_nt1->addItem(
"fptotp",m_fptotp);
382 status = m_nt1->addItem(
"fpt",m_fpt);
383 status = m_nt1->addItem(
"fpte",m_fpte);
384 status = m_nt1->addItem(
"fptmu",m_fptmu);
385 status = m_nt1->addItem(
"fptk",m_fptk);
386 status = m_nt1->addItem(
"fptp",m_fptp);
388 status = m_nt1->addItem(
"lptot",m_lptot);
389 status = m_nt1->addItem(
"lptote",m_lptote);
390 status = m_nt1->addItem(
"lptotmu",m_lptotmu);
391 status = m_nt1->addItem(
"lptotk",m_lptotk);
392 status = m_nt1->addItem(
"lptotp",m_lptotp);
393 status = m_nt1->addItem(
"lpt",m_lpt);
394 status = m_nt1->addItem(
"lpte",m_lpte);
395 status = m_nt1->addItem(
"lptmu",m_lptmu);
396 status = m_nt1->addItem(
"lptk",m_lptk);
397 status = m_nt1->addItem(
"lptp",m_lptp);
399 status = m_nt1->addItem(
"zsigp",m_zsigp);
400 status = m_nt1->addItem(
"zsigpe",m_zsigpe);
401 status = m_nt1->addItem(
"zsigpmu",m_zsigpmu);
402 status = m_nt1->addItem(
"zsigpk",m_zsigpk);
403 status = m_nt1->addItem(
"zsigpp",m_zsigpp);
404 status = m_nt1->addItem(
"fhelix",5,m_fhelix);
405 status = m_nt1->addItem(
"fhelixe",5,m_fhelixe);
406 status = m_nt1->addItem(
"fhelixmu",5,m_fhelixmu);
407 status = m_nt1->addItem(
"fhelixk",5,m_fhelixk);
408 status = m_nt1->addItem(
"fhelixp",5,m_fhelixp);
409 status = m_nt1->addItem(
"lhelix",5,m_lhelix);
410 status = m_nt1->addItem(
"lhelixe",5,m_lhelixe);
411 status = m_nt1->addItem(
"lhelixmu",5,m_lhelixmu);
412 status = m_nt1->addItem(
"lhelixk",5,m_lhelixk);
413 status = m_nt1->addItem(
"lhelixp",5,m_lhelixp);
415 status = m_nt1->addItem(
"zerror",15,m_zerror);
416 status = m_nt1->addItem(
"zerrore",15,m_zerrore);
417 status = m_nt1->addItem(
"zerrormu",15,m_zerrormu);
418 status = m_nt1->addItem(
"zerrork",15,m_zerrork);
419 status = m_nt1->addItem(
"zerrorp",15,m_zerrorp);
420 status = m_nt1->addItem(
"ferror",15,m_ferror);
421 status = m_nt1->addItem(
"ferrore",15,m_ferrore);
422 status = m_nt1->addItem(
"ferrormu",15,m_ferrormu);
423 status = m_nt1->addItem(
"ferrork",15,m_ferrork);
424 status = m_nt1->addItem(
"ferrorp",15,m_ferrorp);
425 status = m_nt1->addItem(
"lerror",15,m_lerror);
426 status = m_nt1->addItem(
"lerrore",15,m_lerrore);
427 status = m_nt1->addItem(
"lerrormu",15,m_lerrormu);
428 status = m_nt1->addItem(
"lerrork",15,m_lerrork);
429 status = m_nt1->addItem(
"lerrorp",15,m_lerrorp);
432 status = m_nt1->addItem(
"evtid",m_evtid);
433 status = m_nt1->addItem(
"mchelix",5,m_mchelix);
434 status = m_nt1->addItem(
"mcptot",m_mcptot);
435 status = m_nt1->addItem(
"mcpid",m_mcpid);
437 if( status.isFailure() ) cout<<
"Ntuple1 add item failed!"<<endl;
443 NTuplePtr nt2(
ntupleSvc(),
"FILE104/n102");
445 if ( nt2 ) m_nt2 = nt2;
447 m_nt2=
ntupleSvc()->book(
"FILE104/n102",CLID_ColumnWiseTuple,
"KalFitComp");
449 status2 = m_nt2->addItem(
"delx",m_delx);
450 status2 = m_nt2->addItem(
"dely",m_dely);
451 status2 = m_nt2->addItem(
"delz",m_delz);
452 status2 = m_nt2->addItem(
"delthe",m_delthe);
453 status2 = m_nt2->addItem(
"delphi",m_delphi);
454 status2 = m_nt2->addItem(
"delp",m_delp);
455 status2 = m_nt2->addItem(
"delpx",m_delpx);
456 status2 = m_nt2->addItem(
"delpy",m_delpy);
457 status2 = m_nt2->addItem(
"delpz",m_delpz);
459 if( status2.isFailure() ) cout<<
"Ntuple2 add item failed!"<<endl;
465 NTuplePtr nt3(
ntupleSvc(),
"FILE104/n103");
467 if ( nt3 ) m_nt3 = nt3;
469 m_nt3=
ntupleSvc()->book(
"FILE104/n103",CLID_ColumnWiseTuple,
"PatRec");
471 status3 = m_nt3->addItem(
"trkhelix",5,m_trkhelix);
472 status3 = m_nt3->addItem(
"trkptot",m_trkptot);
474 status3 = m_nt3->addItem(
"trkerror",15,m_trkerror);
475 status3 = m_nt3->addItem(
"trksigp",m_trksigp);
477 status3 = m_nt3->addItem(
"trkndf",m_trkndf);
478 status3 = m_nt3->addItem(
"trkchisq",m_trkchisq);
479 if( status3.isFailure() ) cout<<
"Ntuple3 add item failed!"<<endl;
485 NTuplePtr nt4(
ntupleSvc(),
"FILE104/n104");
487 if ( nt4 ) m_nt4 = nt4;
489 m_nt4=
ntupleSvc()->book(
"FILE104/n104",CLID_ColumnWiseTuple,
"PatRecComp");
491 status4 = m_nt4->addItem(
"trkdelx",m_trkdelx);
492 status4 = m_nt4->addItem(
"trkdely",m_trkdely);
493 status4 = m_nt4->addItem(
"trkdelz",m_trkdelz);
494 status4 = m_nt4->addItem(
"trkdelthe",m_trkdelthe);
495 status4 = m_nt4->addItem(
"trkdelphi",m_trkdelphi);
496 status4 = m_nt4->addItem(
"trkdelp",m_trkdelp);
497 if( status4.isFailure() ) cout<<
"Ntuple4 add item failed!"<<endl;
502 NTuplePtr nt5(
ntupleSvc(),
"FILE104/n105");
504 if ( nt5 ) m_nt5 = nt5;
506 m_nt5=
ntupleSvc()->book(
"FILE104/n105",CLID_ColumnWiseTuple,
"KalFitdChisq");
508 status5 = m_nt5->addItem(
"dchi2",m_dchi2);
509 status5 = m_nt5->addItem(
"masshyp",m_masshyp);
510 status5 = m_nt5->addItem(
"residual_estim",m_residest);
511 status5 = m_nt5->addItem(
"residual",m_residnew);
512 status5 = m_nt5->addItem(
"layer",m_layer);
513 status5 = m_nt5->addItem(
"kaldr",m_anal_dr);
514 status5 = m_nt5->addItem(
"kalphi0",m_anal_phi0);
515 status5 = m_nt5->addItem(
"kalkappa",m_anal_kappa);
516 status5 = m_nt5->addItem(
"kaldz",m_anal_dz);
517 status5 = m_nt5->addItem(
"kaltanl",m_anal_tanl);
518 status5 = m_nt5->addItem(
"dr_ea",m_anal_ea_dr);
519 status5 = m_nt5->addItem(
"phi0_ea",m_anal_ea_phi0);
520 status5 = m_nt5->addItem(
"kappa_ea",m_anal_ea_kappa);
521 status5 = m_nt5->addItem(
"dz_ea",m_anal_ea_dz);
522 status5 = m_nt5->addItem(
"tanl_ea",m_anal_ea_tanl);
523 if( status5.isFailure() ) cout<<
"Ntuple5 add item failed!"<<endl;
526 NTuplePtr nt6(
ntupleSvc(),
"FILE104/n106");
528 if ( nt6 ) m_nt6 = nt6;
530 m_nt6=
ntupleSvc()->book(
"FILE104/n106",CLID_ColumnWiseTuple,
"kal seg");
532 status6 = m_nt6->addItem(
"docaInc",m_docaInc);
533 status6 = m_nt6->addItem(
"docaExc",m_docaExc);
534 status6 = m_nt6->addItem(
"tdr",m_tdrift);
535 status6 = m_nt6->addItem(
"layerid", m_layerid);
536 status6 = m_nt6->addItem(
"event", m_eventNo);
537 status6 = m_nt6->addItem(
"residualInc",m_residualInc);
538 status6 = m_nt6->addItem(
"residualExc",m_residualExc);
539 status6 = m_nt6->addItem(
"lr",m_lr);
540 status6 = m_nt6->addItem(
"dd",m_dd);
541 status6 = m_nt6->addItem(
"y",m_yposition);
543 if( status6.isFailure() ) cout<<
"Ntuple6 add item failed!"<<endl;
656 MsgStream log(
msgSvc(), name());
657 log << MSG::INFO <<
"in execute()" << endreq;
704 StatusCode sc = service (
"MagneticFieldSvc",IMFSvc);
705 if(sc!=StatusCode::SUCCESS) {
706 log << MSG::ERROR <<
"Unable to open Magnetic field service"<<endreq;
715 std::cout<<
" execute, referField from MagneticFieldSvc: "<< (IMFSvc->
getReferField())*10000 <<std::endl;
722 IDataProviderSvc* evtSvc =
NULL;
723 Gaudi::svcLocator()->service(
"EventDataSvc", evtSvc);
725 log << MSG::INFO <<
"makeTds:event Svc has been found" << endreq;
727 log << MSG::FATAL <<
"makeTds:Could not find eventSvc" << endreq;
728 return StatusCode::SUCCESS;
732 IDataManagerSvc *dataManSvc;
733 dataManSvc=
dynamic_cast<IDataManagerSvc*
>(evtSvc);
734 DataObject *aKalTrackCol;
735 evtSvc->findObject(
"/Event/Recon/RecMdcKalTrackCol",aKalTrackCol);
736 if(aKalTrackCol !=
NULL) {
737 dataManSvc->clearSubTree(
"/Event/Recon/RecMdcKalTrackCol");
738 evtSvc->unregisterObject(
"/Event/Recon/RecMdcKalTrackCol");
741 kalsc = evtSvc->registerObject(
"/Event/Recon/RecMdcKalTrackCol", kalcol);
742 if( kalsc.isFailure() ) {
743 log << MSG::FATAL <<
"Could not register RecMdcKalTrack" << endreq;
744 return StatusCode::SUCCESS;
746 log << MSG::INFO <<
"RecMdcKalTrackCol registered successfully!" <<endreq;
748 DataObject *aKalHelixSegCol;
749 evtSvc->findObject(
"/Event/Recon/RecMdcKalHelixSegCol", aKalHelixSegCol);
750 if(aKalHelixSegCol !=
NULL){
751 dataManSvc->clearSubTree(
"/Event/Recon/RecMdcKalHelixSegCol");
752 evtSvc->unregisterObject(
"/Event/Recon/RecMdcKalHelixSegCol");
755 kalsc = evtSvc->registerObject(
"/Event/Recon/RecMdcKalHelixSegCol", helixsegcol);
756 if( kalsc.isFailure()){
757 log<< MSG::FATAL <<
"Could not register RecMdcKalHelixSeg" <<endreq;
758 return StatusCode::SUCCESS;
760 log << MSG::INFO <<
"RecMdcKalHelixSegCol register successfully!" <<endreq;
772 std::cout<<
"ERROR OCCUR when dynamic_cast in KalFitAlg::execute ...!!"<<std::endl;
775 SmartDataPtr<Event::EventHeader> eventHeader(eventSvc(),
"/Event/EventHeader");
777 log << MSG::WARNING <<
"Could not find Event Header" << endreq;
778 return StatusCode::FAILURE;
780 int eventNo = eventHeader->eventNumber();
781 int runNo = eventHeader->runNumber();
786 SmartDataPtr<RecEsTimeCol> estimeCol(eventSvc(),
"/Event/Recon/RecEsTimeCol");
787 if (estimeCol && estimeCol->size()) {
788 RecEsTimeCol::iterator iter_evt = estimeCol->begin();
789 t0 = (*iter_evt)->getTest();
792 log << MSG::WARNING <<
"Could not find EvTimeCol" << endreq;
793 return StatusCode::SUCCESS;
798 std::cout<<
"in KalFitAlg , we get the event start time = "<<t0<<std::endl;
802 SmartDataPtr<MdcDigiCol> mdcDigiCol(evtSvc,
"/Event/Digi/MdcDigiCol");
803 if (sc!=StatusCode::SUCCESS) {
804 log << MSG::FATAL <<
"Could not find MdcDigiCol!" << endreq;
805 return StatusCode::SUCCESS;
816 m_evtid = eventHeader->eventNumber();
819 SmartDataPtr<McParticleCol> mcPartCol(eventSvc(),
"/Event/MC/McParticleCol");
821 log << MSG::WARNING <<
"Could not find McParticle" << endreq;
826 McParticleCol::iterator i_mcTrk = mcPartCol->begin();
827 for (;i_mcTrk != mcPartCol->end(); i_mcTrk++) {
828 if(!(*i_mcTrk)->primaryParticle())
continue;
829 const HepLorentzVector& mom((*i_mcTrk)->initialFourMomentum());
830 const HepLorentzVector& pos = (*i_mcTrk)->initialPosition();
831 log << MSG::DEBUG <<
"MCINFO:particleId=" << (*i_mcTrk)->particleProperty()
832 <<
" theta=" << mom.theta() <<
" phi="<< mom.phi()
833 <<
" px="<< mom.px() <<
" py="<< mom.py() <<
" pz="<< mom.pz()
836 int pid = (*i_mcTrk)->particleProperty();
838 charge = m_particleTable->particle( pid )->charge();
839 }
else if ( pid <0 ) {
840 charge = m_particleTable->particle( -pid )->charge();
843 log << MSG::WARNING <<
"wrong particle id, please check data" <<endreq;
846 Hep3Vector mom2(mom.px(),mom.py(),mom.pz());
849 log << MSG::DEBUG <<
"charge of the track " <<
charge << endreq;
850 if(
debug_ == 4) cout<<
"helix: "<<mchelix.
a()<<endl;
852 for(
int j =0; j<5; j++) {
853 m_mchelix[j] = mchelix.
a()[j];
856 m_mcptot = sqrt(1+pow(m_mchelix[4],2))/m_mchelix[2];
864 SmartDataPtr<RecMdcTrackCol> newtrkCol(eventSvc(),
"/Event/Recon/RecMdcTrackCol");
866 log << MSG::FATAL <<
"Could not find RecMdcTrackCol" << endreq;
867 return( StatusCode::SUCCESS);
869 log << MSG::INFO <<
"Begin to make MdcRecTrkCol and MdcRecWirhitCol"<<endreq;
874 mtrkadd_mgr->clear();
878 double trkx1,trkx2,trky1,trky2,trkz1,trkz2,trkthe1,trkthe2,trkphi1,trkphi2,trkp1,trkp2,trkr1,trkr2,trkkap1,trkkap2,trktanl1,trktanl2;
882 RecMdcTrackCol::iterator iter_trk = newtrkCol->begin();
883 for(
int kj = 1; iter_trk != newtrkCol->end(); iter_trk++,kj++) {
885 csmp3[kj-1]=(*iter_trk)->p3();
886 csmphi[kj-1] = (*iter_trk)->phi();
890 for(
int j = 0, ij = 0; j<5; j++) {
891 m_trkhelix[j] = (*iter_trk)->helix()[j];
893 for(
int k=0; k<=j; k++,ij++) {
894 m_trkerror[ij] = (*iter_trk)->err()[j][k];
898 m_trkptot = sqrt(1+pow(m_trkhelix[4],2))/m_trkhelix[2];
900 m_trksigp = sqrt(pow((m_trkptot/m_trkhelix[2]),2)*m_trkerror[5]+
901 pow((m_trkhelix[4]/m_trkptot),2)*pow((1/m_trkhelix[2]),4)*m_trkerror[14]-
902 2*m_trkhelix[4]*m_trkerror[12]*pow((1/m_trkhelix[2]),3));
904 m_trkndf = (*iter_trk)->ndof();
905 m_trkchisq = (*iter_trk)->chi2();
907 if (
debug_ == 4) cout<<
"Ea from RecMdcTrackCol..." <<(*iter_trk)->err()<<endl;
909 StatusCode sc3 = m_nt3->write();
910 if( sc3.isFailure() ) cout<<
"Ntuple3 filling failed!"<<endl;
940 log << MSG::DEBUG <<
"retrieved MDC tracks:"
941 <<
" Nhits " <<(*iter_trk)->getNhits()
942 <<
" Nster " <<(*iter_trk)->nster() <<endreq;
944 HitRefVec gothits = (*iter_trk)->getVecHits();
948 rectrk->
id = (*iter_trk)->trackId();
949 rectrk->
chiSq = (*iter_trk)->chi2();
950 rectrk->
ndf = (*iter_trk)->ndof();
951 rectrk->
fiTerm = (*iter_trk)->getFiTerm();
952 rectrk->
nhits = (*iter_trk)->getNhits();
953 rectrk->
nster = (*iter_trk)->nster();
955 rectrk->
stat = (*iter_trk)->stat();
956 status_temp = (*iter_trk)->stat();
958 trkadd->
id = (*iter_trk)->trackId();
962 trkadd->
body = rectrk;
963 rectrk->
add = trkadd;
965 for (
int i=0; i<5; i++) {
966 rectrk->
helix[i] = (*iter_trk)->helix()[i];
967 if( i<3 ) rectrk->
pivot[i] = (*iter_trk)->getPivot()[i];
968 for(
int j = 1; j<i+2;j++) {
969 rectrk->
error[i*(i+1)/2+j-1] = (*iter_trk)->err()(i+1,j);
972 std::sort(gothits.begin(), gothits.end(), order_rechits);
973 HitRefVec::iterator it_gothit = gothits.begin();
974 for( ; it_gothit != gothits.end(); it_gothit++) {
976 if( (*it_gothit)->getStat() != 1 ) {
978 log<<MSG::WARNING<<
"this hit is not used in helix fitting!"<<endreq;
983 log << MSG::DEBUG <<
"retrieved hits in MDC tracks:"
984 <<
" hits DDL " <<(*it_gothit)->getDriftDistLeft()
985 <<
" hits DDR " <<(*it_gothit)->getDriftDistRight()
986 <<
" error DDL " <<(*it_gothit)->getErrDriftDistLeft()
987 <<
" error DDR " <<(*it_gothit)->getErrDriftDistRight()
988 <<
" id of hit "<<(*it_gothit)->getId()
989 <<
" track id of hit "<<(*it_gothit)->getTrkId()
990 <<
" hits ADC " <<(*it_gothit)->getAdc() << endreq;
993 whit->
id = (*it_gothit)->getId();
994 whit->
ddl = (*it_gothit)->getDriftDistLeft();
995 whit->
ddr = (*it_gothit)->getDriftDistRight();
996 whit->
erddl = (*it_gothit)->getErrDriftDistLeft();
997 whit->
erddr = (*it_gothit)->getErrDriftDistRight();
998 whit->
pChiSq = (*it_gothit)->getChisqAdd();
999 whit->
lr = (*it_gothit)->getFlagLR();
1000 whit->
stat = (*it_gothit)->getStat();
1001 mdcid = (*it_gothit)->getMdcId();
1005 int wid = w0id + localwid;
1007 <<
"lr from PR: "<<whit->
lr
1008 <<
" layerId = " << layid
1009 <<
" wireId = " << localwid
1019 whit->
tdc = (*it_gothit)->getTdc();
1020 whit->
adc= (*it_gothit)->getAdc();
1021 rectrk->
hitcol.push_back(whit);
1022 mhit_mgr->push_back(*whit);
1024 mtrk_mgr->push_back(*rectrk);
1025 mtrkadd_mgr->push_back(*trkadd);
1033 m_trkdelx = trkx1 - trkx2;
1034 m_trkdely = trky1 - trky2;
1035 m_trkdelz = trkz1 - trkz2;
1036 m_trkdelthe = trkthe1 + trkthe2;
1037 m_trkdelphi = trkphi1- trkphi2;
1038 m_trkdelp = trkp1 - trkp2;
1039 StatusCode sc4 = m_nt4->write();
1040 if( sc4.isFailure() ) cout<<
"Ntuple4 filling failed!"<<endl;
1043 if(
debug_ == 4) { std::cout<<
"before refit,ntrk,nhits,nadd..."<<mtrk_mgr->size()
1044 <<
"********"<<mhit_mgr->size()<<
"****"<<mtrkadd_mgr->size()<<endl;
1050 double mdang = 180.0 - csmp3[0].angle(csmp3[1].
unit())*180.0/
M_PI;
1051 double mdphi = 180.0 - fabs(csmphi[0]-csmphi[1])*180.0/
M_PI;
1056 log << MSG::DEBUG <<
"after kalman_fitting(),but in execute...."<<endreq;
1063 SmartDataPtr<RecMdcKalTrackCol> recmdckaltrkCol(eventSvc(),
"/Event/Recon/RecMdcKalTrackCol");
1068 RecMdcKalTrackCol::iterator KalTrk= recmdckaltrkCol->begin();
1070 for(;KalTrk !=recmdckaltrkCol->end();KalTrk++){
1072 for(
int hypo=0; hypo<5; hypo++)
1074 if((*KalTrk)->getStat(0,hypo)==1) nFailedTrks[hypo]++;
1077 HelixSegRefVec::iterator iter_hit = gothelixsegs.begin();
1079 int nhitofthistrk=0;
1080 for( ; iter_hit != gothelixsegs.end(); iter_hit++){
1087 iter_hit=gothelixsegs.begin();
1135 return StatusCode::SUCCESS;
3085 MsgStream log(
msgSvc(), name());
3086 double Pt_threshold(0.3);
3087 Hep3Vector IP(0,0,0);
3093 if ( !&whMgr )
return;
3096 int ntrk = mdcMgr->size();
3097 double* rPt =
new double[ntrk];
3098 int* rOM =
new int[ntrk];
3099 unsigned int* order =
new unsigned int[ntrk];
3100 unsigned int* rCont =
new unsigned int[ntrk];
3101 unsigned int* rGen =
new unsigned int[ntrk];
3104 for(vector<MdcRec_trk>::iterator it = mdcMgr->begin(),
3105 end = mdcMgr->end(); it != end; it++) {
3109 rPt[index] = 1 / fabs(it->helix[2]);
3110 if(
debug_ == 4) cout<<
"rPt...[ "<<index<<
" ]...."<< rPt[index] <<endl;
3111 if(rPt[index] < 0) rPt[index] =
DBL_MAX;
3113 std::vector<MdcRec_wirhit*> pt = it->hitcol ;
3114 if(
debug_ == 4) cout<<
"ppt size: "<< pt.size()<<endl;
3116 for (vector<MdcRec_wirhit*>::iterator ii = pt.end()-1;
3117 ii !=pt.begin()-1; ii--) {
3118 int lyr((*ii)->geo->Lyr()->Id());
3119 if (outermost < lyr) outermost = lyr;
3120 if(
debug_ == 4) cout<<
"outmost: "<<outermost<<
" lyr: "<<lyr<<endl;
3122 rOM[index] = outermost;
3123 order[index] = index;
3128 for (
int j, k = ntrk - 1; k >= 0; k = j){
3130 for(
int i = 1; i <= k; i++)
3131 if(rPt[i - 1] < rPt[i]){
3133 std::swap(order[i], order[j]);
3134 std::swap(rPt[i], rPt[j]);
3135 std::swap(rOM[i], rOM[j]);
3136 std::swap(rCont[i], rCont[j]);
3137 std::swap(rGen[i], rGen[j]);
3145 DataObject *aReconEvent;
3146 eventSvc()->findObject(
"/Event/Recon",aReconEvent);
3150 StatusCode sc = eventSvc()->registerObject(
"/Event/Recon",recevt );
3151 if(sc!=StatusCode::SUCCESS) {
3152 log << MSG::FATAL <<
"Could not register ReconEvent" <<endreq;
3160 log << MSG::INFO <<
"beginning to make RecMdcKalTrackCol" <<endreq;
3162 for(
int l = 0; l < ntrk; l++) {
3172 MdcRec_trk& TrasanTRK = *(mdcMgr->begin() + order[l]);
3176 int trasqual = TrasanTRK_add.
quality;
3177 if(
debug_ == 4) cout<<
"kalman_fitting>trasqual... "<<trasqual<<endl;
3178 if (trasqual)
continue;
3182 cout <<
"******* KalFit NUMBER : " << newcount << std::endl;
3186 if ((TrasanTRK_add.
decision & 32) == 32 ||
3187 (TrasanTRK_add.
decision & 64) == 64) type = 1;
3192 TrasanTRK.
pivot[2]);
3195 for(
int i = 0; i < 5; i++)
3196 a[i] = TrasanTRK.
helix[i];
3199 for(
int i = 0, k = 0; i < 5; i++) {
3200 for(
int j = 0; j <= i; j++) {
3202 ea[j][i] = ea[i][j];
3205 double fiTerm = TrasanTRK.
fiTerm;
3213 int inlyr(999), outlyr(-1);
3214 int* rStat =
new int[43];
3215 for(
int irStat=0;irStat<43;++irStat)rStat[irStat]=0;
3216 std::vector<MdcRec_wirhit*> pt=TrasanTRK.
hitcol;
3218 if(
debug_ == 4) cout<<
"*********Pt size****"<< pt.size()<<endl;
3221 for (vector<MdcRec_wirhit*>::iterator ii = pt.end()-1;
3222 ii != pt.begin()-1; ii--) {
3223 Num[(*ii)->geo->Lyr()->Id()]++;
3226 for (vector<MdcRec_wirhit*>::iterator ii = pt.end()-1;
3227 ii != pt.begin()-1; ii--) {
3229 if (Num[(*ii)->geo->Lyr()->Id()]>3) {
3231 cout <<
"WARNING: I found " << Num[(*ii)->geo->Lyr()->Id()]
3232 <<
" hits in the layer "
3233 << (*ii)->geo->Lyr()->Id() << std::endl;
3251 double dist[2] = {rechit.
ddl, rechit.
ddr};
3252 double erdist[2] = {rechit.
erddl, rechit.
erddr};
3257 if (rechit.
lr==2 || rechit.
lr==0) lr_decision=-1;
3259 else if (rechit.
lr==1) lr_decision=1;
3262 int ind(geo->
Lyr()->
Id());
3266 lr_decision, rechit.
tdc,
3271 if (inlyr>ind) inlyr = ind;
3272 if (outlyr<ind) outlyr = ind;
3275 cout <<
"**** NUMBER OF Mdc HITS (TRASAN) = " << hit_asso << std::endl;
3278 int empty_between(0), empty(0);
3279 for (
int i= inlyr; i <= outlyr; i++)
3280 if (!rStat[i]) empty_between++;
3281 empty = empty_between+inlyr+(42-outlyr);
3288 for(std::vector<KalFitHitMdc>::iterator it_hit = track_lead.
HitsMdc().begin(); it_hit!=track_lead.
HitsMdc().end(); it_hit++){
3294 track_lead.
type(type);
3295 unsigned int nhit = track_lead.
HitsMdc().size();
3296 if (!nhit &&
debug_ == 4) {
3297 cout <<
" ATTENTION TRACK WITH ONLY HITS " << nhit << std::endl;
3302 double KalFitst(0), KalFitax(0), KalFitschi2(0);
3305 Hep3Vector outer_pivot(track_lead.
x(fiTerm));
3307 track_lead.
pivot(outer_pivot);
3313 HepSymMatrix Eakal(5,0);
3316 double costheta = track_lead.
a()[4] / sqrt(1.0 + track_lead.
a()[4]*track_lead.
a()[4]);
3326 cout <<
"from Mdc Pattern Recognition: " << std::endl;
3332 cout <<
" dr = " << work.
a()[0]
3333 <<
", Er_dr = " << sqrt(work.
Ea()[0][0]) << std::endl;
3334 cout <<
" phi0 = " << work.
a()[1]
3335 <<
", Er_phi0 = " << sqrt(work.
Ea()[1][1]) << std::endl;
3336 cout <<
" PT = " << 1/work.
a()[2]
3337 <<
", Er_kappa = " << sqrt(work.
Ea()[2][2]) << std::endl;
3338 cout <<
" dz = " << work.
a()[3]
3339 <<
", Er_dz = " << sqrt(work.
Ea()[3][3]) << std::endl;
3340 cout <<
" tanl = " << work.
a()[4]
3341 <<
", Er_tanl = " << sqrt(work.
Ea()[4][4]) << std::endl;
3351 cout <<
" Mdc FIRST KALMAN FIT " << std::endl;
3356 cout <<
" dr = " << work.
a()[0]
3357 <<
", Er_dr = " << sqrt(work.
Ea()[0][0]) << std::endl;
3358 cout <<
" phi0 = " << work.
a()[1]
3359 <<
", Er_phi0 = " << sqrt(work.
Ea()[1][1]) << std::endl;
3360 cout <<
" PT = " << 1/work.
a()[2]
3361 <<
", Er_kappa = " << sqrt(work.
Ea()[2][2]) << std::endl;
3362 cout <<
" dz = " << work.
a()[3]
3363 <<
", Er_dz = " << sqrt(work.
Ea()[3][3]) << std::endl;
3364 cout <<
" tanl = " << work.
a()[4]
3365 <<
", Er_tanl = " << sqrt(work.
Ea()[4][4]) << std::endl;
3371 complete_track(TrasanTRK, TrasanTRK_add, track_lead, kaltrk,kalcol,segcol,1);
3375 IDataProviderSvc* eventSvc =
NULL;
3376 Gaudi::svcLocator()->service(
"EventDataSvc", eventSvc);
3378 log << MSG::INFO <<
"makeTds:event Svc has been found" << endreq;
3380 log << MSG::FATAL <<
"makeTds:Could not find eventSvc" << endreq;
3385 IDataManagerSvc *dataManSvc;
3386 dataManSvc=
dynamic_cast<IDataManagerSvc*
>(eventSvc);
3387 DataObject *aKalTrackCol;
3388 eventSvc->findObject(
"/Event/Recon/RecMdcKalTrackCol",aKalTrackCol);
3389 if(aKalTrackCol !=
NULL) {
3390 dataManSvc->clearSubTree(
"/Event/Recon/RecMdcKalTrackCol");
3391 eventSvc->unregisterObject(
"/Event/Recon/RecMdcKalTrackCol");
3394 kalsc = eventSvc->registerObject(
"/Event/Recon/RecMdcKalTrackCol", kalcol);
3395 if( kalsc.isFailure() ) {
3396 log << MSG::FATAL <<
"Could not register RecMdcKalTrack" << endreq;
3399 log << MSG::INFO <<
"RecMdcKalTrackCol registered successfully!" <<endreq;
3403 DataObject *aRecKalSegEvent;
3404 eventSvc->findObject(
"/Event/Recon/RecMdcKalHelixSegCol", aRecKalSegEvent);
3405 if(aRecKalSegEvent!=
NULL) {
3407 segsc = eventSvc->unregisterObject(
"/Event/Recon/RecMdcKalHelixSegCol");
3408 if(segsc != StatusCode::SUCCESS) {
3409 log << MSG::FATAL <<
"Could not unregister RecMdcKalHelixSegCol collection" << endreq;
3414 segsc = eventSvc->registerObject(
"/Event/Recon/RecMdcKalHelixSegCol", segcol);
3415 if( segsc.isFailure() ) {
3416 log << MSG::FATAL <<
"Could not register RecMdcKalHelixSeg" << endreq;
3419 log << MSG::INFO <<
"RecMdcKalHelixSegCol registered successfully!" <<endreq;
3422 double x1(0.),x2(0.),y1(0.),y2(0.),z1(0.),z2(0.),the1(0.),the2(0.),
phi1(0.),
phi2(0.),
p1(0.),
p2(0.);
3423 double r1(0.),r2(0.),kap1(999.),kap2(999.),tanl1(0),tanl2(0.);
3424 double px1(0.),px2(0.),py1(0.),py2(0.),pz1(0.),pz2(0.),charge1(0.),charge2(0.);
3427 SmartDataPtr<RecMdcKalTrackCol> kaltrkCol(eventSvc,
"/Event/Recon/RecMdcKalTrackCol");
3429 log << MSG::FATAL <<
"Could not find RecMdcKalTrackCol" << endreq;
3432 log << MSG::INFO <<
"Begin to check RecMdcKalTrackCol"<<endreq;
3433 RecMdcKalTrackCol::iterator iter_trk = kaltrkCol->begin();
3434 for(
int jj=1; iter_trk != kaltrkCol->end(); iter_trk++,jj++) {
3435 log << MSG::DEBUG <<
"retrieved MDC Kalmantrack:"
3436 <<
"Track Id: " << (*iter_trk)->getTrackId()
3437 <<
" Mass of the fit: "<< (*iter_trk)->getMass(2)<< endreq
3438 <<
" Length of the track: "<< (*iter_trk)->getLength(2)
3439 <<
" Tof of the track: "<< (*iter_trk)->getTof(2) << endreq
3440 <<
" Chisq of the fit: "<< (*iter_trk)->getChisq(0,2)
3441 <<
" "<< (*iter_trk)->getChisq(1,2) << endreq
3442 <<
"Ndf of the fit: "<< (*iter_trk)->getNdf(0,2)
3443 <<
" "<< (*iter_trk)->getNdf(1,2) << endreq
3444 <<
"Kappa " << (*iter_trk)->getZHelix()[2]
3446 for(
int i = 0; i<43; i++) {
3447 log << MSG::DEBUG <<
"retrieved pathl["<<i<<
"]= "
3448 << (*iter_trk)->getPathl(i) <<endreq;
3452 m_trackid = (*iter_trk)->getTrackId();
3454 for(
int jj =0, iii=0; jj<5; jj++) {
3456 m_length[jj] = (*iter_trk)->getLength(jj);
3457 m_tof[jj] = (*iter_trk)->getTof(jj);
3458 m_nhits[jj] = (*iter_trk)->getNhits(jj);
3459 m_zhelix[jj] = (*iter_trk)->getZHelix()[jj];
3460 m_zhelixe[jj] = (*iter_trk)->getZHelixE()[jj];
3461 m_zhelixmu[jj] = (*iter_trk)->getZHelixMu()[jj];
3462 m_zhelixk[jj] = (*iter_trk)->getZHelixK()[jj];
3463 m_zhelixp[jj] = (*iter_trk)->getZHelixP()[jj];
3464 m_fhelix[jj] = (*iter_trk)->getFHelix()[jj];
3465 m_fhelixe[jj] = (*iter_trk)->getFHelixE()[jj];
3466 m_fhelixmu[jj] = (*iter_trk)->getFHelixMu()[jj];
3467 m_fhelixk[jj] = (*iter_trk)->getFHelixK()[jj];
3468 m_fhelixp[jj] = (*iter_trk)->getFHelixP()[jj];
3469 m_lhelix[jj] = (*iter_trk)->getLHelix()[jj];
3470 m_lhelixe[jj] = (*iter_trk)->getLHelixE()[jj];
3471 m_lhelixmu[jj] = (*iter_trk)->getLHelixMu()[jj];
3472 m_lhelixk[jj] = (*iter_trk)->getLHelixK()[jj];
3473 m_lhelixp[jj] = (*iter_trk)->getLHelixP()[jj];
3476 for(
int kk=0; kk<=jj; kk++,iii++) {
3477 m_zerror[iii] = (*iter_trk)->getZError()[jj][kk];
3478 m_zerrore[iii] = (*iter_trk)->getZErrorE()[jj][kk];
3479 m_zerrormu[iii] = (*iter_trk)->getZErrorMu()[jj][kk];
3480 m_zerrork[iii] = (*iter_trk)->getZErrorK()[jj][kk];
3481 m_zerrorp[iii] = (*iter_trk)->getZErrorP()[jj][kk];
3482 m_ferror[iii] = (*iter_trk)->getFError()[jj][kk];
3483 m_ferrore[iii] = (*iter_trk)->getFErrorE()[jj][kk];
3484 m_ferrormu[iii] = (*iter_trk)->getFErrorMu()[jj][kk];
3485 m_ferrork[iii] = (*iter_trk)->getFErrorK()[jj][kk];
3486 m_ferrorp[iii] = (*iter_trk)->getFErrorP()[jj][kk];
3487 m_lerror[iii] = (*iter_trk)->getLError()[jj][kk];
3488 m_lerrore[iii] = (*iter_trk)->getLErrorE()[jj][kk];
3489 m_lerrormu[iii] = (*iter_trk)->getLErrorMu()[jj][kk];
3490 m_lerrork[iii] = (*iter_trk)->getLErrorK()[jj][kk];
3491 m_lerrorp[iii] = (*iter_trk)->getLErrorP()[jj][kk];
3531 m_chisq[0][0] = (*iter_trk)->getChisq(0,0);
3532 m_chisq[1][0] = (*iter_trk)->getChisq(0,1);
3533 m_chisq[2][0] = (*iter_trk)->getChisq(0,2);
3534 m_chisq[3][0] = (*iter_trk)->getChisq(0,3);
3535 m_chisq[4][0] = (*iter_trk)->getChisq(0,4);
3536 m_chisq[0][1] = (*iter_trk)->getChisq(1,0);
3537 m_chisq[1][1] = (*iter_trk)->getChisq(1,1);
3538 m_chisq[2][1] = (*iter_trk)->getChisq(1,2);
3539 m_chisq[3][1] = (*iter_trk)->getChisq(1,3);
3540 m_chisq[4][1] = (*iter_trk)->getChisq(1,4);
3542 m_ndf[0][0] = (*iter_trk)->getNdf(0,0);
3543 m_ndf[1][0] = (*iter_trk)->getNdf(0,1);
3544 m_ndf[2][0] = (*iter_trk)->getNdf(0,2);
3545 m_ndf[3][0] = (*iter_trk)->getNdf(0,3);
3546 m_ndf[4][0] = (*iter_trk)->getNdf(0,4);
3547 m_ndf[0][1] = (*iter_trk)->getNdf(1,0);
3548 m_ndf[1][1] = (*iter_trk)->getNdf(1,1);
3549 m_ndf[2][1] = (*iter_trk)->getNdf(1,2);
3550 m_ndf[3][1] = (*iter_trk)->getNdf(1,3);
3551 m_ndf[4][1] = (*iter_trk)->getNdf(1,4);
3553 m_stat[0][0] = (*iter_trk)->getStat(0,0);
3554 m_stat[1][0] = (*iter_trk)->getStat(0,1);
3555 m_stat[2][0] = (*iter_trk)->getStat(0,2);
3556 m_stat[3][0] = (*iter_trk)->getStat(0,3);
3557 m_stat[4][0] = (*iter_trk)->getStat(0,4);
3558 m_stat[0][1] = (*iter_trk)->getStat(1,0);
3559 m_stat[1][1] = (*iter_trk)->getStat(1,1);
3560 m_stat[2][1] = (*iter_trk)->getStat(1,2);
3561 m_stat[3][1] = (*iter_trk)->getStat(1,3);
3562 m_stat[4][1] = (*iter_trk)->getStat(1,4);
3564 m_fptot = sqrt(1+pow(m_fhelix[4],2))/m_fhelix[2];
3565 m_fptote = sqrt(1+pow(m_fhelixe[4],2))/m_fhelixe[2];
3566 m_fptotmu = sqrt(1+pow(m_fhelixmu[4],2))/m_fhelixmu[2];
3567 m_fptotk = sqrt(1+pow(m_fhelixk[4],2))/m_fhelixk[2];
3568 m_fptotp = sqrt(1+pow(m_fhelixp[4],2))/m_fhelixp[2];
3570 m_lptot = sqrt(1+pow(m_lhelix[4],2))/m_lhelix[2];
3571 m_lptote = sqrt(1+pow(m_lhelixe[4],2))/m_lhelixe[2];
3572 m_lptotmu = sqrt(1+pow(m_lhelixmu[4],2))/m_lhelixmu[2];
3573 m_lptotk = sqrt(1+pow(m_lhelixk[4],2))/m_lhelixk[2];
3574 m_lptotp = sqrt(1+pow(m_lhelixp[4],2))/m_lhelixp[2];
3576 m_lpt = 1/m_lhelix[2];
3577 m_lpte = 1/m_lhelixe[2];
3578 m_lptmu = 1/m_lhelixmu[2];
3579 m_lptk = 1/m_lhelixk[2];
3580 m_lptp = 1/m_lhelixp[2];
3582 m_fpt = 1/m_fhelix[2];
3583 m_fpte = 1/m_fhelixe[2];
3584 m_fptmu = 1/m_fhelixmu[2];
3585 m_fptk = 1/m_fhelixk[2];
3586 m_fptp = 1/m_fhelixp[2];
3589 std::cout<<
" "<<std::endl;
3590 std::cout<<
"in file Kalman_fitting_anal ,the m_fpt is .."<<m_fpt<<std::endl;
3591 std::cout<<
"in file Kalman_fitting_anal ,the m_fpte is .."<<m_fpte<<std::endl;
3592 std::cout<<
"in file Kalman_fitting_anal ,the m_fptmu is .."<<m_fptmu<<std::endl;
3593 std::cout<<
"in file Kalman_fitting_anal ,the m_fptk is .."<<m_fptk<<std::endl;
3594 std::cout<<
"in file Kalman_fitting_anal ,the m_fptp is .."<<m_fptp<<std::endl;
3597 m_zpt = 1/m_zhelix[2];
3598 m_zpte = 1/m_zhelixe[2];
3599 m_zptmu = 1/m_zhelixmu[2];
3600 m_zptk = 1/m_zhelixk[2];
3601 m_zptp = 1/m_zhelixp[2];
3604 std::cout<<
"in file Kalman_fitting_anal ,the m_zpt is .."<<m_zpt<<std::endl;
3605 std::cout<<
"in file Kalman_fitting_anal ,the m_zpte is .."<<m_zpte<<std::endl;
3606 std::cout<<
"in file Kalman_fitting_anal ,the m_zptmu is .."<<m_zptmu<<std::endl;
3607 std::cout<<
"in file Kalman_fitting_anal ,the m_zptk is .."<<m_zptk<<std::endl;
3608 std::cout<<
"in file Kalman_fitting_anal ,the m_zptp is .."<<m_zptp<<std::endl;
3610 m_zptot = sqrt(1+pow(m_zhelix[4],2))/m_zhelix[2];
3611 m_zptote = sqrt(1+pow(m_zhelixe[4],2))/m_zhelixe[2];
3612 m_zptotmu = sqrt(1+pow(m_zhelixmu[4],2))/m_zhelixmu[2];
3613 m_zptotk = sqrt(1+pow(m_zhelixk[4],2))/m_zhelixk[2];
3614 m_zptotp = sqrt(1+pow(m_zhelixp[4],2))/m_zhelixp[2];
3617 std::cout<<
"in file Kalman_fitting_anal ,the m_zptot is .."<<m_zptot<<std::endl;
3618 std::cout<<
"in file Kalman_fitting_anal ,the m_zptote is .."<<m_zptote<<std::endl;
3619 std::cout<<
"in file Kalman_fitting_anal ,the m_zptotmu is .."<<m_zptotmu<<std::endl;
3620 std::cout<<
"in file Kalman_fitting_anal ,the m_zptotk is .."<<m_zptotk<<std::endl;
3621 std::cout<<
"in file Kalman_fitting_anal ,the m_zptotp is .."<<m_zptotp<<std::endl;
3625 m_zsigp = sqrt(pow((m_zptot/m_zhelix[2]),2)*m_zerror[5]+
3626 pow((m_zhelix[4]/m_zptot),2)*pow((1/m_zhelix[2]),4)*m_zerror[14]-
3627 2*m_zhelix[4]*m_zerror[12]*pow((1/m_zhelix[2]),3));
3628 m_zsigpe = sqrt(pow((m_zptote/m_zhelixe[2]),2)*m_zerrore[5]+
3629 pow((m_zhelixe[4]/m_zptote),2)*pow((1/m_zhelixe[2]),4)*m_zerrore[14]-
3630 2*m_zhelixe[4]*m_zerrore[12]*pow((1/m_zhelixe[2]),3));
3631 m_zsigpmu = sqrt(pow((m_zptotmu/m_zhelixmu[2]),2)*m_zerrormu[5]+
3632 pow((m_zhelixmu[4]/m_zptotmu),2)*pow((1/m_zhelixmu[2]),4)*m_zerrormu[14]-
3633 2*m_zhelixmu[4]*m_zerrormu[12]*pow((1/m_zhelixmu[2]),3));
3634 m_zsigpk = sqrt(pow((m_zptotk/m_zhelixk[2]),2)*m_zerrork[5]+
3635 pow((m_zhelixk[4]/m_zptotk),2)*pow((1/m_zhelixk[2]),4)*m_zerrork[14]-
3636 2*m_zhelixk[4]*m_zerrork[12]*pow((1/m_zhelixk[2]),3));
3637 m_zsigpp = sqrt(pow((m_zptotp/m_zhelixp[2]),2)*m_zerrorp[5]+
3638 pow((m_zhelixp[4]/m_zptotp),2)*pow((1/m_zhelixp[2]),4)*m_zerrorp[14]-
3639 2*m_zhelixp[4]*m_zerrorp[12]*pow((1/m_zhelixp[2]),3));
3642 StatusCode sc1 = m_nt1->write();
3643 if( sc1.isFailure() ) cout<<
"Ntuple1 filling failed!"<<endl;
3648 phi1 = (*iter_trk)->getFFi0();
3649 r1 = (*iter_trk)->getFDr();
3650 z1 = (*iter_trk)->getFDz();
3651 kap1 = (*iter_trk)->getFCpa();
3652 tanl1 = (*iter_trk)->getFTanl();
3653 charge1 = kap1/fabs(kap1);
3656 p1 = sqrt(1+tanl1*tanl1)/kap1;
3657 the1 =
M_PI/2-atan(tanl1);
3660 pz1= tanl1/fabs(kap1);
3662 }
else if(jj == 2) {
3663 phi2 = (*iter_trk)->getFFi0();
3664 r2 = (*iter_trk)->getFDr();
3665 z2 = (*iter_trk)->getFDz();
3666 kap2 = (*iter_trk)->getFCpa();
3667 tanl2 = (*iter_trk)->getFTanl();
3668 charge2 = kap2/fabs(kap2);
3671 p2 = sqrt(1+tanl2*tanl2)/kap1;
3672 the2 =
M_PI/2-atan(tanl2);
3675 pz2= tanl2/fabs(kap2);
3683 m_delthe = the1 + the2;
3686 m_delpx = charge1*fabs(px1) + charge2*fabs(px2);
3687 m_delpy = charge1*fabs(py1) + charge2*fabs(py2);
3688 m_delpz = charge1*fabs(pz1) + charge2*fabs(pz2);
3690 StatusCode sc2 = m_nt2->write();
3691 if( sc2.isFailure() ) cout<<
"Ntuple2 filling failed!"<<endl;
3700 cout <<
"Kalfitting finished " << std::endl;
3706 MsgStream log(
msgSvc(), name());
3707 double Pt_threshold(0.3);
3708 Hep3Vector IP(0,0,0);
3715 if ( !&whMgr )
return;
3718 int ntrk = mdcMgr->size();
3719 double* rPt =
new double[ntrk];
3720 int* rOM =
new int[ntrk];
3721 unsigned int* order =
new unsigned int[ntrk];
3722 unsigned int* rCont =
new unsigned int[ntrk];
3723 unsigned int* rGen =
new unsigned int[ntrk];
3726 for(vector<MdcRec_trk>::iterator it = mdcMgr->begin(),
3727 end = mdcMgr->end(); it != end; it++) {
3731 rPt[index] = 1 / fabs(it->helix[2]);
3732 if(
debug_ == 4) cout<<
"rPt...[ "<<index<<
" ]...."<< rPt[index] <<endl;
3733 if(rPt[index] < 0) rPt[index] =
DBL_MAX;
3735 std::vector<MdcRec_wirhit*> pt = it->hitcol ;
3736 if(
debug_ == 4) cout<<
"ppt size: "<< pt.size()<<endl;
3738 for (vector<MdcRec_wirhit*>::iterator ii = pt.end()-1;
3739 ii !=pt.begin()-1; ii--) {
3740 int lyr((*ii)->geo->Lyr()->Id());
3741 if (outermost < lyr) outermost = lyr;
3742 if(
debug_ == 4) cout<<
"outmost: "<<outermost<<
" lyr: "<<lyr<<endl;
3744 rOM[index] = outermost;
3745 order[index] = index;
3750 for (
int j, k = ntrk - 1; k >= 0; k = j){
3752 for(
int i = 1; i <= k; i++)
3753 if(rPt[i - 1] < rPt[i]){
3755 std::swap(order[i], order[j]);
3756 std::swap(rPt[i], rPt[j]);
3757 std::swap(rOM[i], rOM[j]);
3758 std::swap(rCont[i], rCont[j]);
3759 std::swap(rGen[i], rGen[j]);
3766 DataObject *aReconEvent;
3767 eventSvc()->findObject(
"/Event/Recon",aReconEvent);
3771 StatusCode sc = eventSvc()->registerObject(
"/Event/Recon",recevt );
3772 if(sc!=StatusCode::SUCCESS) {
3773 log << MSG::FATAL <<
"Could not register ReconEvent" <<endreq;
3781 log << MSG::INFO <<
"beginning to make RecMdcKalTrackCol" <<endreq;
3784 for(
int l = 0; l < ntrk; l++) {
3786 MdcRec_trk& TrasanTRK = *(mdcMgr->begin() + order[l]);
3790 int trasqual = TrasanTRK_add.
quality;
3791 if(
debug_ == 4) cout<<
"kalman_fitting>trasqual... "<<trasqual<<endl;
3792 if (trasqual)
continue;
3796 cout <<
"******* KalFit NUMBER : " << newcount << std::endl;
3800 if ((TrasanTRK_add.
decision & 32) == 32 ||
3801 (TrasanTRK_add.
decision & 64) == 64) type = 1;
3806 TrasanTRK.
pivot[2]);
3809 for(
int i = 0; i < 5; i++)
3810 a[i] = TrasanTRK.
helix[i];
3813 for(
int i = 0, k = 0; i < 5; i++) {
3814 for(
int j = 0; j <= i; j++) {
3816 ea[j][i] = ea[i][j];
3822 double fiTerm = TrasanTRK.
fiTerm;
3828 int inlyr(999), outlyr(-1);
3829 int* rStat =
new int[43];
3830 for(
int irStat=0;irStat<43;++irStat)rStat[irStat]=0;
3831 std::vector<MdcRec_wirhit*> pt=TrasanTRK.
hitcol;
3833 if(
debug_ == 4) cout<<
"*********Pt size****"<< pt.size()<<endl;
3836 for (vector<MdcRec_wirhit*>::iterator ii = pt.end()-1;
3837 ii != pt.begin()-1; ii--) {
3838 Num[(*ii)->geo->Lyr()->Id()]++;
3842 for (vector<MdcRec_wirhit*>::iterator ii = pt.end()-1;
3843 ii != pt.begin()-1; ii--) {
3846 if (Num[(*ii)->geo->Lyr()->Id()]>3) {
3848 cout <<
"WARNING: I found " << Num[(*ii)->geo->Lyr()->Id()]
3849 <<
" hits in the layer "
3850 << (*ii)->geo->Lyr()->Id() << std::endl;
3877 double dist[2] = {rechit.
ddl, rechit.
ddr};
3878 double erdist[2] = {rechit.
erddl, rechit.
erddr};
3883 if (rechit.
lr==2 || rechit.
lr==0) lr_decision=-1;
3885 else if (rechit.
lr==1) lr_decision=1;
3888 int ind(geo->
Lyr()->
Id());
3890 lr_decision, rechit.
tdc,
3895 if (inlyr>ind) inlyr = ind;
3896 if (outlyr<ind) outlyr = ind;
3899 cout <<
"**** NUMBER OF Mdc HITS (TRASAN) = " << hit_asso << std::endl;
3902 int empty_between(0), empty(0);
3903 for (
int i= inlyr; i <= outlyr; i++)
3904 if (!rStat[i]) empty_between++;
3905 empty = empty_between+inlyr+(42-outlyr);
3910 track_lead.
type(type);
3911 unsigned int nhit = track_lead.
HitsMdc().size();
3912 if (!nhit &&
debug_ == 4) {
3913 cout <<
" ATTENTION TRACK WITH ONLY HITS " << nhit << std::endl;
3918 double KalFitst(0), KalFitax(0), KalFitschi2(0);
3920 Hep3Vector outer_pivot(track_lead.
x(fiTerm));
3923 std::cout<<
"before track_lead.pivot(outer_pivot) ,the error matrix of track_lead is .."<<track_lead.
Ea()<<std::endl;
3925 track_lead.
pivot(outer_pivot);
3930 HepSymMatrix Eakal(5,0);
3934 double costheta = track_lead.
a()[4] / sqrt(1.0 + track_lead.
a()[4]*track_lead.
a()[4]);
3944 std::cout <<
"from Mdc Pattern Recognition: " << std::endl;
3950 std::cout <<
" dr = " << work.
a()[0]
3951 <<
", Er_dr = " << sqrt(work.
Ea()[0][0]) << std::endl;
3952 std::cout <<
" phi0 = " << work.
a()[1]
3953 <<
", Er_phi0 = " << sqrt(work.
Ea()[1][1]) << std::endl;
3954 std::cout <<
" PT = " << 1/work.
a()[2]
3955 <<
", Er_kappa = " << sqrt(work.
Ea()[2][2]) << std::endl;
3956 std::cout <<
" dz = " << work.
a()[3]
3957 <<
", Er_dz = " << sqrt(work.
Ea()[3][3]) << std::endl;
3958 std::cout <<
" tanl = " << work.
a()[4]
3959 <<
", Er_tanl = " << sqrt(work.
Ea()[4][4]) << std::endl;
3967 cout <<
" Mdc FIRST KALMAN FIT " << std::endl;
3972 cout <<
" dr = " << work.
a()[0]
3973 <<
", Er_dr = " << sqrt(work.
Ea()[0][0]) << std::endl;
3974 cout <<
" phi0 = " << work.
a()[1]
3975 <<
", Er_phi0 = " << sqrt(work.
Ea()[1][1]) << std::endl;
3976 cout <<
" PT = " << 1/work.
a()[2]
3977 <<
", Er_kappa = " << sqrt(work.
Ea()[2][2]) << std::endl;
3978 cout <<
" dz = " << work.
a()[3]
3979 <<
", Er_dz = " << sqrt(work.
Ea()[3][3]) << std::endl;
3980 cout <<
" tanl = " << work.
a()[4]
3981 <<
", Er_tanl = " << sqrt(work.
Ea()[4][4]) << std::endl;
3988 complete_track(TrasanTRK, TrasanTRK_add, track_lead, kaltrk,kalcol,segcol);
3994 DataObject *aRecKalEvent;
3995 eventSvc()->findObject(
"/Event/Recon/RecMdcKalTrackCol", aRecKalEvent);
3996 if(aRecKalEvent!=
NULL) {
3998 kalsc = eventSvc()->unregisterObject(
"/Event/Recon/RecMdcKalTrackCol");
3999 if(kalsc != StatusCode::SUCCESS) {
4000 log << MSG::FATAL <<
"Could not unregister RecMdcKalTrack collection" << endreq;
4005 kalsc = eventSvc()->registerObject(
"/Event/Recon/RecMdcKalTrackCol", kalcol);
4006 if( kalsc.isFailure()) {
4007 log << MSG::FATAL <<
"Could not register RecMdcKalTrack" << endreq;
4010 log << MSG::INFO <<
"RecMdcKalTrackCol registered successfully!" <<endreq;
4016 DataObject *aRecKalSegEvent;
4017 eventSvc()->findObject(
"/Event/Recon/RecMdcKalHelixSegCol", aRecKalSegEvent);
4018 if(aRecKalSegEvent!=
NULL) {
4020 segsc = eventSvc()->unregisterObject(
"/Event/Recon/RecMdcKalHelixSegCol");
4021 if(segsc != StatusCode::SUCCESS) {
4022 log << MSG::FATAL <<
"Could not unregister RecMdcKalHelixSegCol collection" << endreq;
4027 segsc = eventSvc()->registerObject(
"/Event/Recon/RecMdcKalHelixSegCol", segcol);
4028 if( segsc.isFailure() ) {
4029 log << MSG::FATAL <<
"Could not register RecMdcKalHelixSeg" << endreq;
4032 log << MSG::INFO <<
"RecMdcKalHelixSegCol registered successfully!" <<endreq;
4035 double x1(0.),x2(0.),y1(0.),y2(0.),z1(0.),z2(0.),the1(0.),the2(0.),
phi1(0.),
phi2(0.),
p1(0.),
p2(0.);
4036 double r1(0.),r2(0.),kap1(999.),kap2(999.),tanl1(0.),tanl2(0.);
4039 SmartDataPtr<RecMdcKalTrackCol> kaltrkCol(eventSvc(),
"/Event/Recon/RecMdcKalTrackCol");
4041 log << MSG::FATAL <<
"Could not find RecMdcKalTrackCol" << endreq;
4044 log << MSG::INFO <<
"Begin to check RecMdcKalTrackCol"<<endreq;
4045 RecMdcKalTrackCol::iterator iter_trk = kaltrkCol->begin();
4046 for(
int jj=1; iter_trk != kaltrkCol->end(); iter_trk++,jj++) {
4047 log << MSG::DEBUG <<
"retrieved MDC Kalmantrack:"
4048 <<
"Track Id: " << (*iter_trk)->getTrackId()
4049 <<
" Mass of the fit: "<< (*iter_trk)->getMass(2)<< endreq
4050 <<
" Length of the track: "<< (*iter_trk)->getLength(2)
4051 <<
" Tof of the track: "<< (*iter_trk)->getTof(2) << endreq
4052 <<
" Chisq of the fit: "<< (*iter_trk)->getChisq(0,2)
4053 <<
" "<< (*iter_trk)->getChisq(1,2) << endreq
4054 <<
"Ndf of the fit: "<< (*iter_trk)->getNdf(0,1)
4055 <<
" "<< (*iter_trk)->getNdf(1,2) << endreq
4056 <<
"Kappa " << (*iter_trk)->getZHelix()[2]
4061 std::cout<<
"the size of gothelixsegs ..."<<gothelixsegs.size()<<std::endl;
4064 HelixSegRefVec::iterator it_gothelixseg = gothelixsegs.begin();
4065 for( ; it_gothelixseg != gothelixsegs.end(); it_gothelixseg++) {
4067 std::cout<<
"the layerId of this helixseg is ..."<<(*it_gothelixseg)->getLayerId()<<std::endl;
4068 std::cout<<
"the residual of this helixseg exclude the meas hit"<<(*it_gothelixseg)->getResExcl()<<std::endl;
4069 std::cout<<
"the residual of this helixseg include the meas hit"<<(*it_gothelixseg)->getResIncl()<<std::endl;
4070 std::cout<<
"the track id of the helixseg is ..."<<(*it_gothelixseg)->getTrackId() <<std::endl;
4071 std::cout<<
"the tof of the helixseg is ..."<<(*it_gothelixseg)->getTof()<<std::endl;
4072 std::cout<<
"the Zhit of the helixseg is ..."<<(*it_gothelixseg)->getZhit()<<std::endl;
4075 for(
int i = 0; i<43; i++) {
4076 log << MSG::DEBUG <<
"retrieved pathl["<<i<<
"]= "
4077 << (*iter_trk)->getPathl(i) <<endreq;
4081 m_trackid = (*iter_trk)->getTrackId();
4082 for(
int jj =0, iii=0; jj<5; jj++) {
4083 m_length[jj] = (*iter_trk)->getLength(jj);
4084 m_tof[jj] = (*iter_trk)->getTof(jj);
4085 m_nhits[jj] = (*iter_trk)->getNhits(jj);
4086 m_zhelix[jj] = (*iter_trk)->getZHelix()[jj];
4087 m_zhelixe[jj] = (*iter_trk)->getZHelixE()[jj];
4088 m_zhelixmu[jj] = (*iter_trk)->getZHelixMu()[jj];
4089 m_zhelixk[jj] = (*iter_trk)->getZHelixK()[jj];
4090 m_zhelixp[jj] = (*iter_trk)->getZHelixP()[jj];
4091 m_fhelix[jj] = (*iter_trk)->getFHelix()[jj];
4092 m_fhelixe[jj] = (*iter_trk)->getFHelixE()[jj];
4093 m_fhelixmu[jj] = (*iter_trk)->getFHelixMu()[jj];
4094 m_fhelixk[jj] = (*iter_trk)->getFHelixK()[jj];
4095 m_fhelixp[jj] = (*iter_trk)->getFHelixP()[jj];
4096 m_lhelix[jj] = (*iter_trk)->getLHelix()[jj];
4097 m_lhelixe[jj] = (*iter_trk)->getLHelixE()[jj];
4098 m_lhelixmu[jj] = (*iter_trk)->getLHelixMu()[jj];
4099 m_lhelixk[jj] = (*iter_trk)->getLHelixK()[jj];
4100 m_lhelixp[jj] = (*iter_trk)->getLHelixP()[jj];
4102 for(
int kk=0; kk<=jj; kk++,iii++) {
4103 m_zerror[iii] = (*iter_trk)->getZError()[jj][kk];
4104 m_zerrore[iii] = (*iter_trk)->getZErrorE()[jj][kk];
4105 m_zerrormu[iii] = (*iter_trk)->getZErrorMu()[jj][kk];
4106 m_zerrork[iii] = (*iter_trk)->getZErrorK()[jj][kk];
4107 m_zerrorp[iii] = (*iter_trk)->getZErrorP()[jj][kk];
4108 m_ferror[iii] = (*iter_trk)->getFError()[jj][kk];
4109 m_ferrore[iii] = (*iter_trk)->getFErrorE()[jj][kk];
4110 m_ferrormu[iii] = (*iter_trk)->getFErrorMu()[jj][kk];
4111 m_ferrork[iii] = (*iter_trk)->getFErrorK()[jj][kk];
4112 m_ferrorp[iii] = (*iter_trk)->getFErrorP()[jj][kk];
4113 m_lerror[iii] = (*iter_trk)->getLError()[jj][kk];
4114 m_lerrore[iii] = (*iter_trk)->getLErrorE()[jj][kk];
4115 m_lerrormu[iii] = (*iter_trk)->getLErrorMu()[jj][kk];
4116 m_lerrork[iii] = (*iter_trk)->getLErrorK()[jj][kk];
4117 m_lerrorp[iii] = (*iter_trk)->getLErrorP()[jj][kk];
4158 m_chisq[0][0] = (*iter_trk)->getChisq(0,0);
4159 m_chisq[1][0] = (*iter_trk)->getChisq(0,1);
4160 m_chisq[2][0] = (*iter_trk)->getChisq(0,2);
4161 m_chisq[3][0] = (*iter_trk)->getChisq(0,3);
4162 m_chisq[4][0] = (*iter_trk)->getChisq(0,4);
4163 m_chisq[0][1] = (*iter_trk)->getChisq(1,0);
4164 m_chisq[1][1] = (*iter_trk)->getChisq(1,1);
4165 m_chisq[2][1] = (*iter_trk)->getChisq(1,2);
4166 m_chisq[3][1] = (*iter_trk)->getChisq(1,3);
4167 m_chisq[4][1] = (*iter_trk)->getChisq(1,4);
4169 m_ndf[0][0] = (*iter_trk)->getNdf(0,0);
4170 m_ndf[1][0] = (*iter_trk)->getNdf(0,1);
4171 m_ndf[2][0] = (*iter_trk)->getNdf(0,2);
4172 m_ndf[3][0] = (*iter_trk)->getNdf(0,3);
4173 m_ndf[4][0] = (*iter_trk)->getNdf(0,4);
4174 m_ndf[0][1] = (*iter_trk)->getNdf(1,0);
4175 m_ndf[1][1] = (*iter_trk)->getNdf(1,1);
4176 m_ndf[2][1] = (*iter_trk)->getNdf(1,2);
4177 m_ndf[3][1] = (*iter_trk)->getNdf(1,3);
4178 m_ndf[4][1] = (*iter_trk)->getNdf(1,4);
4180 m_stat[0][0] = (*iter_trk)->getStat(0,0);
4181 m_stat[1][0] = (*iter_trk)->getStat(0,1);
4182 m_stat[2][0] = (*iter_trk)->getStat(0,2);
4183 m_stat[3][0] = (*iter_trk)->getStat(0,3);
4184 m_stat[4][0] = (*iter_trk)->getStat(0,4);
4185 m_stat[0][1] = (*iter_trk)->getStat(1,0);
4186 m_stat[1][1] = (*iter_trk)->getStat(1,1);
4187 m_stat[2][1] = (*iter_trk)->getStat(1,2);
4188 m_stat[3][1] = (*iter_trk)->getStat(1,3);
4189 m_stat[4][1] = (*iter_trk)->getStat(1,4);
4191 m_fptot = sqrt(1+pow(m_fhelix[4],2))/m_fhelix[2];
4192 m_fptote = sqrt(1+pow(m_fhelixe[4],2))/m_fhelixe[2];
4193 m_fptotmu = sqrt(1+pow(m_fhelixmu[4],2))/m_fhelixmu[2];
4194 m_fptotk = sqrt(1+pow(m_fhelixk[4],2))/m_fhelixk[2];
4195 m_fptotp = sqrt(1+pow(m_fhelixp[4],2))/m_fhelixp[2];
4197 m_zpt = 1/m_zhelix[2];
4198 m_zpte = 1/m_zhelixe[2];
4199 m_zptmu = 1/m_zhelixmu[2];
4200 m_zptk = 1/m_zhelixk[2];
4201 m_zptp = 1/m_zhelixp[2];
4203 m_fpt = 1/m_fhelix[2];
4204 m_fpte = 1/m_fhelixe[2];
4205 m_fptmu = 1/m_fhelixmu[2];
4206 m_fptk = 1/m_fhelixk[2];
4207 m_fptp = 1/m_fhelixp[2];
4209 m_lpt = 1/m_lhelix[2];
4210 m_lpte = 1/m_lhelixe[2];
4211 m_lptmu = 1/m_lhelixmu[2];
4212 m_lptk = 1/m_lhelixk[2];
4213 m_lptp = 1/m_lhelixp[2];
4215 m_lptot = sqrt(1+pow(m_lhelix[4],2))/m_lhelix[2];
4216 m_lptote = sqrt(1+pow(m_lhelixe[4],2))/m_lhelixe[2];
4217 m_lptotmu = sqrt(1+pow(m_lhelixmu[4],2))/m_lhelixmu[2];
4218 m_lptotk = sqrt(1+pow(m_lhelixk[4],2))/m_lhelixk[2];
4219 m_lptotp = sqrt(1+pow(m_lhelixp[4],2))/m_lhelixp[2];
4221 m_zptot = sqrt(1+pow(m_zhelix[4],2))/m_zhelix[2];
4222 m_zptote = sqrt(1+pow(m_zhelixe[4],2))/m_zhelixe[2];
4223 m_zptotmu = sqrt(1+pow(m_zhelixmu[4],2))/m_zhelixmu[2];
4224 m_zptotk = sqrt(1+pow(m_zhelixk[4],2))/m_zhelixk[2];
4225 m_zptotp = sqrt(1+pow(m_zhelixp[4],2))/m_zhelixp[2];
4227 m_zsigp = sqrt(pow((m_zptot/m_zhelix[2]),2)*m_zerror[5]+
4228 pow((m_zhelix[4]/m_zptot),2)*pow((1/m_zhelix[2]),4)*m_zerror[14]-
4229 2*m_zhelix[4]*m_zerror[12]*pow((1/m_zhelix[2]),3));
4230 m_zsigpe = sqrt(pow((m_zptote/m_zhelixe[2]),2)*m_zerrore[5]+
4231 pow((m_zhelixe[4]/m_zptote),2)*pow((1/m_zhelixe[2]),4)*m_zerrore[14]-
4232 2*m_zhelixe[4]*m_zerrore[12]*pow((1/m_zhelixe[2]),3));
4233 m_zsigpmu = sqrt(pow((m_zptotmu/m_zhelixmu[2]),2)*m_zerrormu[5]+
4234 pow((m_zhelixmu[4]/m_zptotmu),2)*pow((1/m_zhelixmu[2]),4)*m_zerrormu[14]-
4235 2*m_zhelixmu[4]*m_zerrormu[12]*pow((1/m_zhelixmu[2]),3));
4236 m_zsigpk = sqrt(pow((m_zptotk/m_zhelixk[2]),2)*m_zerrork[5]+
4237 pow((m_zhelixk[4]/m_zptotk),2)*pow((1/m_zhelixk[2]),4)*m_zerrork[14]-
4238 2*m_zhelixk[4]*m_zerrork[12]*pow((1/m_zhelixk[2]),3));
4239 m_zsigpp = sqrt(pow((m_zptotp/m_zhelixp[2]),2)*m_zerrorp[5]+
4240 pow((m_zhelixp[4]/m_zptotp),2)*pow((1/m_zhelixp[2]),4)*m_zerrorp[14]-
4241 2*m_zhelixp[4]*m_zerrorp[12]*pow((1/m_zhelixp[2]),3));
4244 StatusCode sc1 = m_nt1->write();
4245 if( sc1.isFailure() ) cout<<
"Ntuple1 filling failed!"<<endl;
4250 phi1 = (*iter_trk)->getFFi0();
4251 r1 = (*iter_trk)->getFDr();
4252 z1 = (*iter_trk)->getFDz();
4253 kap1 = (*iter_trk)->getFCpa();
4254 tanl1 = (*iter_trk)->getFTanl();
4257 p1 = sqrt(1+tanl1*tanl1)/kap1;
4258 the1 =
M_PI/2-atan(tanl1);
4259 }
else if(jj == 2) {
4260 phi2 = (*iter_trk)->getFFi0();
4261 r2 = (*iter_trk)->getFDr();
4262 z2 = (*iter_trk)->getFDz();
4263 kap2 = (*iter_trk)->getFCpa();
4264 tanl2 = (*iter_trk)->getFTanl();
4267 p2 = sqrt(1+tanl2*tanl2)/kap1;
4268 the2 =
M_PI/2-atan(tanl2);
4276 m_delthe = the1 + the2;
4279 StatusCode sc2 = m_nt2->write();
4280 if( sc2.isFailure() ) cout<<
"Ntuple2 filling failed!"<<endl;
4288 cout <<
"Kalfitting finished " << std::endl;
4295 MsgStream log(
msgSvc(), name());
4296 double Pt_threshold(0.3);
4297 Hep3Vector IP(0,0,0);
4304 if ( !&whMgr )
return;
4307 int ntrk = mdcMgr->size();
4310 int nhits = whMgr->size();
4315 DataObject *aReconEvent;
4316 eventSvc()->findObject(
"/Event/Recon",aReconEvent);
4320 StatusCode sc = eventSvc()->registerObject(
"/Event/Recon",recevt );
4321 if(sc!=StatusCode::SUCCESS) {
4322 log << MSG::FATAL <<
"Could not register ReconEvent" <<endreq;
4330 log << MSG::INFO <<
"beginning to make RecMdcKalTrackCol" <<endreq;
4342 if ((TrasanTRK_add.
decision & 32) == 32 ||
4343 (TrasanTRK_add.
decision & 64) == 64) type = 1;
4348 TrasanTRK.
pivot[2]);
4351 for(
int i = 0; i < 5; i++)
4352 a[i] = TrasanTRK.
helix[i];
4355 for(
int i = 0, k = 0; i < 5; i++) {
4356 for(
int j = 0; j <= i; j++) {
4358 ea[j][i] = ea[i][j];
4364 double fiTerm = TrasanTRK.
fiTerm;
4372 int trasqual = TrasanTRK_add.
quality;
4373 if(
debug_ == 4) cout<<
"kalman_fitting>trasqual... "<<trasqual<<endl;
4374 if (trasqual)
return;
4376 int inlyr(999), outlyr(-1);
4377 int* rStat =
new int[43];
4378 for(
int irStat=0;irStat<43;++irStat)rStat[irStat]=0;
4379 std::vector<MdcRec_wirhit*> pt=TrasanTRK.
hitcol;
4381 if(
debug_ == 4) cout<<
"*********Pt size****"<< pt.size()<<endl;
4384 for (vector<MdcRec_wirhit*>::iterator ii = pt.end()-1;
4385 ii != pt.begin()-1; ii--) {
4386 Num[(*ii)->geo->Lyr()->Id()]++;
4389 for (vector<MdcRec_wirhit*>::iterator ii = pt.end()-1;
4390 ii != pt.begin()-1; ii--) {
4393 if (Num[(*ii)->geo->Lyr()->Id()]>3) {
4395 cout <<
"WARNING: I found " << Num[(*ii)->geo->Lyr()->Id()]
4396 <<
" hits in the layer "
4397 << (*ii)->geo->Lyr()->Id() << std::endl;
4403 double dist[2] = {rechit.
ddl, rechit.
ddr};
4404 double erdist[2] = {rechit.
erddl, rechit.
erddr};
4409 if (rechit.
lr==2 || rechit.
lr==0) lr_decision=-1;
4411 else if (rechit.
lr==1) lr_decision=1;
4414 int ind(geo->
Lyr()->
Id());
4416 lr_decision, rechit.
tdc,
4421 if (inlyr>ind) inlyr = ind;
4422 if (outlyr<ind) outlyr = ind;
4425 int empty_between(0), empty(0);
4426 for (
int i= inlyr; i <= outlyr; i++)
4427 if (!rStat[i]) empty_between++;
4428 empty = empty_between+inlyr+(42-outlyr);
4432 cout <<
"**** NUMBER OF Mdc HITS (TRASAN) = " << hit_asso << std::endl;
4437 track_lead.
type(type);
4438 unsigned int nhit = track_lead.
HitsMdc().size();
4440 cout <<
" ATTENTION TRACK WITH ONLY HITS " << nhit << std::endl;
4445 double KalFitst(0), KalFitax(0), KalFitschi2(0);
4447 Hep3Vector outer_pivot(track_lead.
x(fiTerm));
4450 std::cout<<
"before track_lead.pivot(outer_pivot) ,the error matrix of track_lead is .."<<track_lead.
Ea()<<std::endl;
4452 track_lead.
pivot(outer_pivot);
4457 HepSymMatrix Eakal(5,0);
4461 double costheta = track_lead.
a()[4] / sqrt(1.0 + track_lead.
a()[4]*track_lead.
a()[4]);
4471 std::cout <<
"from Mdc Pattern Recognition: " << std::endl;
4477 std::cout <<
" dr = " << work.
a()[0]
4478 <<
", Er_dr = " << sqrt(work.
Ea()[0][0]) << std::endl;
4479 std::cout <<
" phi0 = " << work.
a()[1]
4480 <<
", Er_phi0 = " << sqrt(work.
Ea()[1][1]) << std::endl;
4481 std::cout <<
" PT = " << 1/work.
a()[2]
4482 <<
", Er_kappa = " << sqrt(work.
Ea()[2][2]) << std::endl;
4483 std::cout <<
" dz = " << work.
a()[3]
4484 <<
", Er_dz = " << sqrt(work.
Ea()[3][3]) << std::endl;
4485 std::cout <<
" tanl = " << work.
a()[4]
4486 <<
", Er_tanl = " << sqrt(work.
Ea()[4][4]) << std::endl;
4493 cout <<
" Mdc FIRST KALMAN FIT " << std::endl;
4498 cout <<
" dr = " << work1.
a()[0]
4499 <<
", Er_dr = " << sqrt(work1.
Ea()[0][0]) << std::endl;
4500 cout <<
" phi0 = " << work1.
a()[1]
4501 <<
", Er_phi0 = " << sqrt(work1.
Ea()[1][1]) << std::endl;
4502 cout <<
" PT = " << 1/work1.
a()[2]
4503 <<
", Er_kappa = " << sqrt(work1.
Ea()[2][2]) << std::endl;
4504 cout <<
" dz = " << work1.
a()[3]
4505 <<
", Er_dz = " << sqrt(work1.
Ea()[3][3]) << std::endl;
4506 cout <<
" tanl = " << work1.
a()[4]
4507 <<
", Er_tanl = " << sqrt(work1.
Ea()[4][4]) << std::endl;
4514 complete_track(TrasanTRK, TrasanTRK_add, track_lead, kaltrk,kalcol,segcol);
4519 DataObject *aRecKalEvent;
4520 eventSvc()->findObject(
"/Event/Recon/RecMdcKalTrackCol", aRecKalEvent);
4521 if(aRecKalEvent!=
NULL) {
4523 kalsc = eventSvc()->unregisterObject(
"/Event/Recon/RecMdcKalTrackCol");
4524 if(kalsc != StatusCode::SUCCESS) {
4525 log << MSG::FATAL <<
"Could not unregister RecMdcKalTrack collection" << endreq;
4530 kalsc = eventSvc()->registerObject(
"/Event/Recon/RecMdcKalTrackCol", kalcol);
4531 if( kalsc.isFailure()) {
4532 log << MSG::FATAL <<
"Could not register RecMdcKalTrack" << endreq;
4535 log << MSG::INFO <<
"RecMdcKalTrackCol registered successfully!" <<endreq;
4541 DataObject *aRecKalSegEvent;
4542 eventSvc()->findObject(
"/Event/Recon/RecMdcKalHelixSegCol", aRecKalSegEvent);
4543 if(aRecKalSegEvent!=
NULL) {
4545 segsc = eventSvc()->unregisterObject(
"/Event/Recon/RecMdcKalHelixSegCol");
4546 if(segsc != StatusCode::SUCCESS) {
4547 log << MSG::FATAL <<
"Could not unregister RecMdcKalHelixSegCol collection" << endreq;
4552 segsc = eventSvc()->registerObject(
"/Event/Recon/RecMdcKalHelixSegCol", segcol);
4553 if( segsc.isFailure() ) {
4554 log << MSG::FATAL <<
"Could not register RecMdcKalHelixSeg" << endreq;
4557 log << MSG::INFO <<
"RecMdcKalHelixSegCol registered successfully!" <<endreq;
4560 double x1(0.),x2(0.),y1(0.),y2(0.),z1(0.),z2(0.),the1(0.),the2(0.),
phi1(0.),
phi2(0.),
p1(0.),
p2(0.);
4561 double r1(0.),r2(0.),kap1(999.),kap2(999.),tanl1(0.),tanl2(0.);
4564 SmartDataPtr<RecMdcKalTrackCol> kaltrkCol(eventSvc(),
"/Event/Recon/RecMdcKalTrackCol");
4566 log << MSG::FATAL <<
"Could not find RecMdcKalTrackCol" << endreq;
4569 log << MSG::INFO <<
"Begin to check RecMdcKalTrackCol"<<endreq;
4570 RecMdcKalTrackCol::iterator iter_trk = kaltrkCol->begin();
4571 for(
int jj=1; iter_trk != kaltrkCol->end(); iter_trk++,jj++) {
4572 log << MSG::DEBUG <<
"retrieved MDC Kalmantrack:"
4573 <<
"Track Id: " << (*iter_trk)->getTrackId()
4574 <<
" Mass of the fit: "<< (*iter_trk)->getMass(2)<< endreq
4575 <<
" Length of the track: "<< (*iter_trk)->getLength(2)
4576 <<
" Tof of the track: "<< (*iter_trk)->getTof(2) << endreq
4577 <<
" Chisq of the fit: "<< (*iter_trk)->getChisq(0,2)
4578 <<
" "<< (*iter_trk)->getChisq(1,2) << endreq
4579 <<
"Ndf of the fit: "<< (*iter_trk)->getNdf(0,1)
4580 <<
" "<< (*iter_trk)->getNdf(1,2) << endreq
4581 <<
"Kappa " << (*iter_trk)->getZHelix()[2]
4582 <<
"zhelixmu "<<(*iter_trk)->getZHelixMu()
4587 std::cout<<
"the size of gothelixsegs ..."<<gothelixsegs.size()<<std::endl;
4590 HelixSegRefVec::iterator it_gothelixseg = gothelixsegs.begin();
4591 for( ; it_gothelixseg != gothelixsegs.end(); it_gothelixseg++) {
4593 std::cout<<
"the layerId of this helixseg is ..."<<(*it_gothelixseg)->getLayerId()<<std::endl;
4594 std::cout<<
"the residual of this helixseg exclude the meas hit"<<(*it_gothelixseg)->getResExcl()<<std::endl;
4595 std::cout<<
"the residual of this helixseg include the meas hit"<<(*it_gothelixseg)->getResIncl()<<std::endl;
4596 std::cout<<
"the track id of the helixseg is ..."<<(*it_gothelixseg)->getTrackId() <<std::endl;
4597 std::cout<<
"the tof of the helixseg is ..."<<(*it_gothelixseg)->getTof()<<std::endl;
4598 std::cout<<
"the Zhit of the helixseg is ..."<<(*it_gothelixseg)->getZhit()<<std::endl;
4601 for(
int i = 0; i<43; i++) {
4602 log << MSG::DEBUG <<
"retrieved pathl["<<i<<
"]= "
4603 << (*iter_trk)->getPathl(i) <<endreq;
4607 m_trackid = (*iter_trk)->getTrackId();
4608 for(
int jj =0, iii=0; jj<5; jj++) {
4609 m_length[jj] = (*iter_trk)->getLength(jj);
4610 m_tof[jj] = (*iter_trk)->getTof(jj);
4611 m_nhits[jj] = (*iter_trk)->getNhits(jj);
4612 m_zhelix[jj] = (*iter_trk)->getZHelix()[jj];
4613 m_zhelixe[jj] = (*iter_trk)->getZHelixE()[jj];
4614 m_zhelixmu[jj] = (*iter_trk)->getZHelixMu()[jj];
4615 m_zhelixk[jj] = (*iter_trk)->getZHelixK()[jj];
4616 m_zhelixp[jj] = (*iter_trk)->getZHelixP()[jj];
4617 m_fhelix[jj] = (*iter_trk)->getFHelix()[jj];
4618 m_fhelixe[jj] = (*iter_trk)->getFHelixE()[jj];
4619 m_fhelixmu[jj] = (*iter_trk)->getFHelixMu()[jj];
4620 m_fhelixk[jj] = (*iter_trk)->getFHelixK()[jj];
4621 m_fhelixp[jj] = (*iter_trk)->getFHelixP()[jj];
4622 m_lhelix[jj] = (*iter_trk)->getLHelix()[jj];
4623 m_lhelixe[jj] = (*iter_trk)->getLHelixE()[jj];
4624 m_lhelixmu[jj] = (*iter_trk)->getLHelixMu()[jj];
4625 m_lhelixk[jj] = (*iter_trk)->getLHelixK()[jj];
4626 m_lhelixp[jj] = (*iter_trk)->getLHelixP()[jj];
4628 for(
int kk=0; kk<=jj; kk++,iii++) {
4629 m_zerror[iii] = (*iter_trk)->getZError()[jj][kk];
4630 m_zerrore[iii] = (*iter_trk)->getZErrorE()[jj][kk];
4631 m_zerrormu[iii] = (*iter_trk)->getZErrorMu()[jj][kk];
4632 m_zerrork[iii] = (*iter_trk)->getZErrorK()[jj][kk];
4633 m_zerrorp[iii] = (*iter_trk)->getZErrorP()[jj][kk];
4634 m_ferror[iii] = (*iter_trk)->getFError()[jj][kk];
4635 m_ferrore[iii] = (*iter_trk)->getFErrorE()[jj][kk];
4636 m_ferrormu[iii] = (*iter_trk)->getFErrorMu()[jj][kk];
4637 m_ferrork[iii] = (*iter_trk)->getFErrorK()[jj][kk];
4638 m_ferrorp[iii] = (*iter_trk)->getFErrorP()[jj][kk];
4639 m_lerror[iii] = (*iter_trk)->getLError()[jj][kk];
4640 m_lerrore[iii] = (*iter_trk)->getLErrorE()[jj][kk];
4641 m_lerrormu[iii] = (*iter_trk)->getLErrorMu()[jj][kk];
4642 m_lerrork[iii] = (*iter_trk)->getLErrorK()[jj][kk];
4643 m_lerrorp[iii] = (*iter_trk)->getLErrorP()[jj][kk];
4650 m_chisq[0][0] = (*iter_trk)->getChisq(0,0);
4651 m_chisq[1][0] = (*iter_trk)->getChisq(0,1);
4652 m_chisq[2][0] = (*iter_trk)->getChisq(0,2);
4653 m_chisq[3][0] = (*iter_trk)->getChisq(0,3);
4654 m_chisq[4][0] = (*iter_trk)->getChisq(0,4);
4655 m_chisq[0][1] = (*iter_trk)->getChisq(1,0);
4656 m_chisq[1][1] = (*iter_trk)->getChisq(1,1);
4657 m_chisq[2][1] = (*iter_trk)->getChisq(1,2);
4658 m_chisq[3][1] = (*iter_trk)->getChisq(1,3);
4659 m_chisq[4][1] = (*iter_trk)->getChisq(1,4);
4661 m_ndf[0][0] = (*iter_trk)->getNdf(0,0);
4662 m_ndf[1][0] = (*iter_trk)->getNdf(0,1);
4663 m_ndf[2][0] = (*iter_trk)->getNdf(0,2);
4664 m_ndf[3][0] = (*iter_trk)->getNdf(0,3);
4665 m_ndf[4][0] = (*iter_trk)->getNdf(0,4);
4666 m_ndf[0][1] = (*iter_trk)->getNdf(1,0);
4667 m_ndf[1][1] = (*iter_trk)->getNdf(1,1);
4668 m_ndf[2][1] = (*iter_trk)->getNdf(1,2);
4669 m_ndf[3][1] = (*iter_trk)->getNdf(1,3);
4670 m_ndf[4][1] = (*iter_trk)->getNdf(1,4);
4672 m_stat[0][0] = (*iter_trk)->getStat(0,0);
4673 m_stat[1][0] = (*iter_trk)->getStat(0,1);
4674 m_stat[2][0] = (*iter_trk)->getStat(0,2);
4675 m_stat[3][0] = (*iter_trk)->getStat(0,3);
4676 m_stat[4][0] = (*iter_trk)->getStat(0,4);
4677 m_stat[0][1] = (*iter_trk)->getStat(1,0);
4678 m_stat[1][1] = (*iter_trk)->getStat(1,1);
4679 m_stat[2][1] = (*iter_trk)->getStat(1,2);
4680 m_stat[3][1] = (*iter_trk)->getStat(1,3);
4681 m_stat[4][1] = (*iter_trk)->getStat(1,4);
4683 m_fptot = sqrt(1+pow(m_fhelix[4],2))/m_fhelix[2];
4684 m_fptote = sqrt(1+pow(m_fhelixe[4],2))/m_fhelixe[2];
4685 m_fptotmu = sqrt(1+pow(m_fhelixmu[4],2))/m_fhelixmu[2];
4686 m_fptotk = sqrt(1+pow(m_fhelixk[4],2))/m_fhelixk[2];
4687 m_fptotp = sqrt(1+pow(m_fhelixp[4],2))/m_fhelixp[2];
4689 m_zpt = 1/m_zhelix[2];
4690 m_zpte = 1/m_zhelixe[2];
4691 m_zptmu = 1/m_zhelixmu[2];
4692 m_zptk = 1/m_zhelixk[2];
4693 m_zptp = 1/m_zhelixp[2];
4695 m_fpt = 1/m_fhelix[2];
4696 m_fpte = 1/m_fhelixe[2];
4697 m_fptmu = 1/m_fhelixmu[2];
4698 m_fptk = 1/m_fhelixk[2];
4699 m_fptp = 1/m_fhelixp[2];
4701 m_lpt = 1/m_lhelix[2];
4702 m_lpte = 1/m_lhelixe[2];
4703 m_lptmu = 1/m_lhelixmu[2];
4704 m_lptk = 1/m_lhelixk[2];
4705 m_lptp = 1/m_lhelixp[2];
4707 m_lptot = sqrt(1+pow(m_lhelix[4],2))/m_lhelix[2];
4708 m_lptote = sqrt(1+pow(m_lhelixe[4],2))/m_lhelixe[2];
4709 m_lptotmu = sqrt(1+pow(m_lhelixmu[4],2))/m_lhelixmu[2];
4710 m_lptotk = sqrt(1+pow(m_lhelixk[4],2))/m_lhelixk[2];
4711 m_lptotp = sqrt(1+pow(m_lhelixp[4],2))/m_lhelixp[2];
4713 m_zptot = sqrt(1+pow(m_zhelix[4],2))/m_zhelix[2];
4714 m_zptote = sqrt(1+pow(m_zhelixe[4],2))/m_zhelixe[2];
4715 m_zptotmu = sqrt(1+pow(m_zhelixmu[4],2))/m_zhelixmu[2];
4716 m_zptotk = sqrt(1+pow(m_zhelixk[4],2))/m_zhelixk[2];
4717 m_zptotp = sqrt(1+pow(m_zhelixp[4],2))/m_zhelixp[2];
4719 m_zsigp = sqrt(pow((m_zptot/m_zhelix[2]),2)*m_zerror[5]+
4720 pow((m_zhelix[4]/m_zptot),2)*pow((1/m_zhelix[2]),4)*m_zerror[14]-
4721 2*m_zhelix[4]*m_zerror[12]*pow((1/m_zhelix[2]),3));
4722 m_zsigpe = sqrt(pow((m_zptote/m_zhelixe[2]),2)*m_zerrore[5]+
4723 pow((m_zhelixe[4]/m_zptote),2)*pow((1/m_zhelixe[2]),4)*m_zerrore[14]-
4724 2*m_zhelixe[4]*m_zerrore[12]*pow((1/m_zhelixe[2]),3));
4725 m_zsigpmu = sqrt(pow((m_zptotmu/m_zhelixmu[2]),2)*m_zerrormu[5]+
4726 pow((m_zhelixmu[4]/m_zptotmu),2)*pow((1/m_zhelixmu[2]),4)*m_zerrormu[14]-
4727 2*m_zhelixmu[4]*m_zerrormu[12]*pow((1/m_zhelixmu[2]),3));
4728 m_zsigpk = sqrt(pow((m_zptotk/m_zhelixk[2]),2)*m_zerrork[5]+
4729 pow((m_zhelixk[4]/m_zptotk),2)*pow((1/m_zhelixk[2]),4)*m_zerrork[14]-
4730 2*m_zhelixk[4]*m_zerrork[12]*pow((1/m_zhelixk[2]),3));
4731 m_zsigpp = sqrt(pow((m_zptotp/m_zhelixp[2]),2)*m_zerrorp[5]+
4732 pow((m_zhelixp[4]/m_zptotp),2)*pow((1/m_zhelixp[2]),4)*m_zerrorp[14]-
4733 2*m_zhelixp[4]*m_zerrorp[12]*pow((1/m_zhelixp[2]),3));
4736 StatusCode sc1 = m_nt1->write();
4737 if( sc1.isFailure() ) cout<<
"Ntuple1 filling failed!"<<endl;
4742 phi1 = (*iter_trk)->getFFi0();
4743 r1 = (*iter_trk)->getFDr();
4744 z1 = (*iter_trk)->getFDz();
4745 kap1 = (*iter_trk)->getFCpa();
4746 tanl1 = (*iter_trk)->getFTanl();
4749 p1 = sqrt(1+tanl1*tanl1)/kap1;
4750 the1 =
M_PI/2-atan(tanl1);
4751 }
else if(jj == 2) {
4752 phi2 = (*iter_trk)->getFFi0();
4753 r2 = (*iter_trk)->getFDr();
4754 z2 = (*iter_trk)->getFDz();
4755 kap2 = (*iter_trk)->getFCpa();
4756 tanl2 = (*iter_trk)->getFTanl();
4759 p2 = sqrt(1+tanl2*tanl2)/kap1;
4760 the2 =
M_PI/2-atan(tanl2);
4768 m_delthe = the1 + the2;
4771 StatusCode sc2 = m_nt2->write();
4772 if( sc2.isFailure() ) cout<<
"Ntuple2 filling failed!"<<endl;
4776 cout <<
"Kalfitting finished " << std::endl;
4784 MsgStream log(
msgSvc(), name());
4785 double Pt_threshold(0.3);
4786 Hep3Vector IP(0,0,0);
4793 if ( !&whMgr )
return;
4796 int ntrk = mdcMgr->size();
4799 int nhits = whMgr->size();
4803 double* rY =
new double[ntrk];
4804 double* rfiTerm =
new double[ntrk];
4805 double* rPt =
new double[ntrk];
4806 int* rOM =
new int[ntrk];
4807 unsigned int* order =
new unsigned int[ntrk];
4808 unsigned int* rCont =
new unsigned int[ntrk];
4809 unsigned int* rGen =
new unsigned int[ntrk];
4812 Hep3Vector csmp3[2];
4813 for(vector<MdcRec_trk>::iterator it = mdcMgr->begin(),
4814 end = mdcMgr->end(); it != end; it++) {
4816 rfiTerm[index]=it->fiTerm;
4821 rPt[index] = 1 / fabs(it->helix[2]);
4822 if(
debug_ == 4) cout<<
"rPt...[ "<<index<<
" ]...."<< rPt[index] <<endl;
4823 if(rPt[index] < 0) rPt[index] =
DBL_MAX;
4825 std::vector<MdcRec_wirhit*> pt = it->hitcol ;
4826 if(
debug_ == 4) cout<<
"ppt size: "<< pt.size()<<endl;
4828 for (vector<MdcRec_wirhit*>::iterator ii = pt.end()-1;
4829 ii !=pt.begin()-1; ii--) {
4830 int lyr((*ii)->geo->Lyr()->Id());
4831 if (outermost < lyr) {
4833 rY[index] = (*ii)->geo->Forward().y();
4835 if(
debug_ == 4) cout<<
"outmost: "<<outermost<<
" lyr: "<<lyr<<endl;
4837 rOM[index] = outermost;
4838 order[index] = index;
4843 for (
int j, k = ntrk - 1; k >= 0; k = j){
4845 for(
int i = 1; i <= k; i++)
4846 if(rY[i - 1] < rY[i]){
4848 std::swap(order[i], order[j]);
4849 std::swap(rY[i], rY[j]);
4850 std::swap(rOM[i], rOM[j]);
4851 std::swap(rCont[i], rCont[j]);
4852 std::swap(rGen[i], rGen[j]);
4861 DataObject *aReconEvent;
4862 eventSvc()->findObject(
"/Event/Recon",aReconEvent);
4866 StatusCode sc = eventSvc()->registerObject(
"/Event/Recon",recevt );
4867 if(sc!=StatusCode::SUCCESS) {
4868 log << MSG::FATAL <<
"Could not register ReconEvent" <<endreq;
4876 log << MSG::INFO <<
"beginning to make RecMdcKalTrackCol" <<endreq;
4883 MdcRec_trk& TrasanTRK = *(mdcMgr->begin() + order[1]);
4892 cout <<
"******* KalFit NUMBER : " << newcount << std::endl;
4896 if ((TrasanTRK_add.
decision & 32) == 32 ||
4897 (TrasanTRK_add.
decision & 64) == 64) type = 1;
4902 TrasanTRK.
pivot[2]);
4905 for(
int i = 0; i < 5; i++)
4906 a[i] = TrasanTRK.
helix[i];
4909 for(
int i = 0, k = 0; i < 5; i++) {
4910 for(
int j = 0; j <= i; j++) {
4912 ea[j][i] = ea[i][j];
4918 double fiTerm = TrasanTRK.
fiTerm;
4925 for(
int l = 0; l < ntrk; l++) {
4926 MdcRec_trk& TrasanTRK1 = *(mdcMgr->begin() + order[l]);
4927 MdcRec_trk_add& TrasanTRK_add1 = *(mdc_addMgr->begin()+order[l]);
4929 int trasqual = TrasanTRK_add1.
quality;
4930 if(
debug_ == 4) cout<<
"kalman_fitting>trasqual... "<<trasqual<<endl;
4931 if (trasqual)
continue;
4933 int inlyr(999), outlyr(-1);
4934 int* rStat =
new int[43];
4935 for(
int irStat=0;irStat<43;++irStat)rStat[irStat]=0;
4936 std::vector<MdcRec_wirhit*> pt=TrasanTRK1.
hitcol;
4938 if(
debug_ == 4) cout<<
"*********Pt size****"<< pt.size()<<endl;
4941 for (vector<MdcRec_wirhit*>::iterator ii = pt.end()-1;
4942 ii != pt.begin()-1; ii--) {
4943 Num[(*ii)->geo->Lyr()->Id()]++;
4946 for (vector<MdcRec_wirhit*>::iterator ii = pt.end()-1;
4947 ii != pt.begin()-1; ii--) {
4950 if (Num[(*ii)->geo->Lyr()->Id()]>3) {
4952 cout <<
"WARNING: I found " << Num[(*ii)->geo->Lyr()->Id()]
4953 <<
" hits in the layer "
4954 << (*ii)->geo->Lyr()->Id() << std::endl;
4960 double dist[2] = {rechit.
ddl, rechit.
ddr};
4961 double erdist[2] = {rechit.
erddl, rechit.
erddr};
4966 if (rechit.
lr==2 || rechit.
lr==0) lr_decision=-1;
4968 else if (rechit.
lr==1) lr_decision=1;
4971 int ind(geo->
Lyr()->
Id());
4973 lr_decision, rechit.
tdc,
4978 if (inlyr>ind) inlyr = ind;
4979 if (outlyr<ind) outlyr = ind;
4982 int empty_between(0), empty(0);
4983 for (
int i= inlyr; i <= outlyr; i++)
4984 if (!rStat[i]) empty_between++;
4985 empty = empty_between+inlyr+(42-outlyr);
4989 cout <<
"**** NUMBER OF Mdc HITS (TRASAN) = " << hit_asso << std::endl;
4994 track_lead.
type(type);
4995 unsigned int nhit = track_lead.
HitsMdc().size();
4997 cout <<
" ATTENTION TRACK WITH ONLY HITS " << nhit << std::endl;
5002 double KalFitst(0), KalFitax(0), KalFitschi2(0);
5004 Hep3Vector outer_pivot(track_lead.
x(fiTerm));
5007 std::cout<<
"before track_lead.pivot(outer_pivot) ,the error matrix of track_lead is .."<<track_lead.
Ea()<<std::endl;
5009 track_lead.
pivot(outer_pivot);
5014 HepSymMatrix Eakal(5,0);
5018 double costheta = track_lead.
a()[4] / sqrt(1.0 + track_lead.
a()[4]*track_lead.
a()[4]);
5028 std::cout <<
"from Mdc Pattern Recognition: " << std::endl;
5034 std::cout <<
" dr = " << work.
a()[0]
5035 <<
", Er_dr = " << sqrt(work.
Ea()[0][0]) << std::endl;
5036 std::cout <<
" phi0 = " << work.
a()[1]
5037 <<
", Er_phi0 = " << sqrt(work.
Ea()[1][1]) << std::endl;
5038 std::cout <<
" PT = " << 1/work.
a()[2]
5039 <<
", Er_kappa = " << sqrt(work.
Ea()[2][2]) << std::endl;
5040 std::cout <<
" dz = " << work.
a()[3]
5041 <<
", Er_dz = " << sqrt(work.
Ea()[3][3]) << std::endl;
5042 std::cout <<
" tanl = " << work.
a()[4]
5043 <<
", Er_tanl = " << sqrt(work.
Ea()[4][4]) << std::endl;
5050 cout <<
" Mdc FIRST KALMAN FIT " << std::endl;
5055 cout <<
" dr = " << work1.
a()[0]
5056 <<
", Er_dr = " << sqrt(work1.
Ea()[0][0]) << std::endl;
5057 cout <<
" phi0 = " << work1.
a()[1]
5058 <<
", Er_phi0 = " << sqrt(work1.
Ea()[1][1]) << std::endl;
5059 cout <<
" PT = " << 1/work1.
a()[2]
5060 <<
", Er_kappa = " << sqrt(work1.
Ea()[2][2]) << std::endl;
5061 cout <<
" dz = " << work1.
a()[3]
5062 <<
", Er_dz = " << sqrt(work1.
Ea()[3][3]) << std::endl;
5063 cout <<
" tanl = " << work1.
a()[4]
5064 <<
", Er_tanl = " << sqrt(work1.
Ea()[4][4]) << std::endl;
5071 complete_track(TrasanTRK, TrasanTRK_add, track_lead, kaltrk,kalcol,segcol);
5077 DataObject *aRecKalEvent;
5078 eventSvc()->findObject(
"/Event/Recon/RecMdcKalTrackCol", aRecKalEvent);
5079 if(aRecKalEvent!=
NULL) {
5081 kalsc = eventSvc()->unregisterObject(
"/Event/Recon/RecMdcKalTrackCol");
5082 if(kalsc != StatusCode::SUCCESS) {
5083 log << MSG::FATAL <<
"Could not unregister RecMdcKalTrack collection" << endreq;
5088 kalsc = eventSvc()->registerObject(
"/Event/Recon/RecMdcKalTrackCol", kalcol);
5089 if( kalsc.isFailure()) {
5090 log << MSG::FATAL <<
"Could not register RecMdcKalTrack" << endreq;
5093 log << MSG::INFO <<
"RecMdcKalTrackCol registered successfully!" <<endreq;
5099 DataObject *aRecKalSegEvent;
5100 eventSvc()->findObject(
"/Event/Recon/RecMdcKalHelixSegCol", aRecKalSegEvent);
5101 if(aRecKalSegEvent!=
NULL) {
5103 segsc = eventSvc()->unregisterObject(
"/Event/Recon/RecMdcKalHelixSegCol");
5104 if(segsc != StatusCode::SUCCESS) {
5105 log << MSG::FATAL <<
"Could not unregister RecMdcKalHelixSegCol collection" << endreq;
5110 segsc = eventSvc()->registerObject(
"/Event/Recon/RecMdcKalHelixSegCol", segcol);
5111 if( segsc.isFailure() ) {
5112 log << MSG::FATAL <<
"Could not register RecMdcKalHelixSeg" << endreq;
5115 log << MSG::INFO <<
"RecMdcKalHelixSegCol registered successfully!" <<endreq;
5118 double x1(0.),x2(0.),y1(0.),y2(0.),z1(0.),z2(0.),the1(0.),the2(0.),
phi1(0.),
phi2(0.),
p1(0.),
p2(0.);
5119 double r1(0.),r2(0.),kap1(999.),kap2(999.),tanl1(0.),tanl2(0.);
5122 SmartDataPtr<RecMdcKalTrackCol> kaltrkCol(eventSvc(),
"/Event/Recon/RecMdcKalTrackCol");
5124 log << MSG::FATAL <<
"Could not find RecMdcKalTrackCol" << endreq;
5127 log << MSG::INFO <<
"Begin to check RecMdcKalTrackCol"<<endreq;
5128 RecMdcKalTrackCol::iterator iter_trk = kaltrkCol->begin();
5129 for(
int jj=1; iter_trk != kaltrkCol->end(); iter_trk++,jj++) {
5130 log << MSG::DEBUG <<
"retrieved MDC Kalmantrack:"
5131 <<
"Track Id: " << (*iter_trk)->getTrackId()
5132 <<
" Mass of the fit: "<< (*iter_trk)->getMass(2)<< endreq
5133 <<
" Length of the track: "<< (*iter_trk)->getLength(2)
5134 <<
" Tof of the track: "<< (*iter_trk)->getTof(2) << endreq
5135 <<
" Chisq of the fit: "<< (*iter_trk)->getChisq(0,2)
5136 <<
" "<< (*iter_trk)->getChisq(1,2) << endreq
5137 <<
"Ndf of the fit: "<< (*iter_trk)->getNdf(0,1)
5138 <<
" "<< (*iter_trk)->getNdf(1,2) << endreq
5139 <<
"Kappa " << (*iter_trk)->getZHelix()[2]
5140 <<
"zhelixmu "<<(*iter_trk)->getZHelixMu()
5145 std::cout<<
"the size of gothelixsegs ..."<<gothelixsegs.size()<<std::endl;
5148 HelixSegRefVec::iterator it_gothelixseg = gothelixsegs.begin();
5149 for( ; it_gothelixseg != gothelixsegs.end(); it_gothelixseg++) {
5151 std::cout<<
"the layerId of this helixseg is ..."<<(*it_gothelixseg)->getLayerId()<<std::endl;
5152 std::cout<<
"the residual of this helixseg exclude the meas hit"<<(*it_gothelixseg)->getResExcl()<<std::endl;
5153 std::cout<<
"the residual of this helixseg include the meas hit"<<(*it_gothelixseg)->getResIncl()<<std::endl;
5154 std::cout<<
"the track id of the helixseg is ..."<<(*it_gothelixseg)->getTrackId() <<std::endl;
5155 std::cout<<
"the tof of the helixseg is ..."<<(*it_gothelixseg)->getTof()<<std::endl;
5156 std::cout<<
"the Zhit of the helixseg is ..."<<(*it_gothelixseg)->getZhit()<<std::endl;
5159 for(
int i = 0; i<43; i++) {
5160 log << MSG::DEBUG <<
"retrieved pathl["<<i<<
"]= "
5161 << (*iter_trk)->getPathl(i) <<endreq;
5165 m_trackid = (*iter_trk)->getTrackId();
5166 for(
int jj =0, iii=0; jj<5; jj++) {
5167 m_length[jj] = (*iter_trk)->getLength(jj);
5168 m_tof[jj] = (*iter_trk)->getTof(jj);
5169 m_nhits[jj] = (*iter_trk)->getNhits(jj);
5170 m_zhelix[jj] = (*iter_trk)->getZHelix()[jj];
5171 m_zhelixe[jj] = (*iter_trk)->getZHelixE()[jj];
5172 m_zhelixmu[jj] = (*iter_trk)->getZHelixMu()[jj];
5173 m_zhelixk[jj] = (*iter_trk)->getZHelixK()[jj];
5174 m_zhelixp[jj] = (*iter_trk)->getZHelixP()[jj];
5175 m_fhelix[jj] = (*iter_trk)->getFHelix()[jj];
5176 m_fhelixe[jj] = (*iter_trk)->getFHelixE()[jj];
5177 m_fhelixmu[jj] = (*iter_trk)->getFHelixMu()[jj];
5178 m_fhelixk[jj] = (*iter_trk)->getFHelixK()[jj];
5179 m_fhelixp[jj] = (*iter_trk)->getFHelixP()[jj];
5180 m_lhelix[jj] = (*iter_trk)->getLHelix()[jj];
5181 m_lhelixe[jj] = (*iter_trk)->getLHelixE()[jj];
5182 m_lhelixmu[jj] = (*iter_trk)->getLHelixMu()[jj];
5183 m_lhelixk[jj] = (*iter_trk)->getLHelixK()[jj];
5184 m_lhelixp[jj] = (*iter_trk)->getLHelixP()[jj];
5186 for(
int kk=0; kk<=jj; kk++,iii++) {
5187 m_zerror[iii] = (*iter_trk)->getZError()[jj][kk];
5188 m_zerrore[iii] = (*iter_trk)->getZErrorE()[jj][kk];
5189 m_zerrormu[iii] = (*iter_trk)->getZErrorMu()[jj][kk];
5190 m_zerrork[iii] = (*iter_trk)->getZErrorK()[jj][kk];
5191 m_zerrorp[iii] = (*iter_trk)->getZErrorP()[jj][kk];
5192 m_ferror[iii] = (*iter_trk)->getFError()[jj][kk];
5193 m_ferrore[iii] = (*iter_trk)->getFErrorE()[jj][kk];
5194 m_ferrormu[iii] = (*iter_trk)->getFErrorMu()[jj][kk];
5195 m_ferrork[iii] = (*iter_trk)->getFErrorK()[jj][kk];
5196 m_ferrorp[iii] = (*iter_trk)->getFErrorP()[jj][kk];
5197 m_lerror[iii] = (*iter_trk)->getLError()[jj][kk];
5198 m_lerrore[iii] = (*iter_trk)->getLErrorE()[jj][kk];
5199 m_lerrormu[iii] = (*iter_trk)->getLErrorMu()[jj][kk];
5200 m_lerrork[iii] = (*iter_trk)->getLErrorK()[jj][kk];
5201 m_lerrorp[iii] = (*iter_trk)->getLErrorP()[jj][kk];
5208 m_chisq[0][0] = (*iter_trk)->getChisq(0,0);
5209 m_chisq[1][0] = (*iter_trk)->getChisq(0,1);
5210 m_chisq[2][0] = (*iter_trk)->getChisq(0,2);
5211 m_chisq[3][0] = (*iter_trk)->getChisq(0,3);
5212 m_chisq[4][0] = (*iter_trk)->getChisq(0,4);
5213 m_chisq[0][1] = (*iter_trk)->getChisq(1,0);
5214 m_chisq[1][1] = (*iter_trk)->getChisq(1,1);
5215 m_chisq[2][1] = (*iter_trk)->getChisq(1,2);
5216 m_chisq[3][1] = (*iter_trk)->getChisq(1,3);
5217 m_chisq[4][1] = (*iter_trk)->getChisq(1,4);
5219 m_ndf[0][0] = (*iter_trk)->getNdf(0,0);
5220 m_ndf[1][0] = (*iter_trk)->getNdf(0,1);
5221 m_ndf[2][0] = (*iter_trk)->getNdf(0,2);
5222 m_ndf[3][0] = (*iter_trk)->getNdf(0,3);
5223 m_ndf[4][0] = (*iter_trk)->getNdf(0,4);
5224 m_ndf[0][1] = (*iter_trk)->getNdf(1,0);
5225 m_ndf[1][1] = (*iter_trk)->getNdf(1,1);
5226 m_ndf[2][1] = (*iter_trk)->getNdf(1,2);
5227 m_ndf[3][1] = (*iter_trk)->getNdf(1,3);
5228 m_ndf[4][1] = (*iter_trk)->getNdf(1,4);
5230 m_stat[0][0] = (*iter_trk)->getStat(0,0);
5231 m_stat[1][0] = (*iter_trk)->getStat(0,1);
5232 m_stat[2][0] = (*iter_trk)->getStat(0,2);
5233 m_stat[3][0] = (*iter_trk)->getStat(0,3);
5234 m_stat[4][0] = (*iter_trk)->getStat(0,4);
5235 m_stat[0][1] = (*iter_trk)->getStat(1,0);
5236 m_stat[1][1] = (*iter_trk)->getStat(1,1);
5237 m_stat[2][1] = (*iter_trk)->getStat(1,2);
5238 m_stat[3][1] = (*iter_trk)->getStat(1,3);
5239 m_stat[4][1] = (*iter_trk)->getStat(1,4);
5241 m_fptot = sqrt(1+pow(m_fhelix[4],2))/m_fhelix[2];
5242 m_fptote = sqrt(1+pow(m_fhelixe[4],2))/m_fhelixe[2];
5243 m_fptotmu = sqrt(1+pow(m_fhelixmu[4],2))/m_fhelixmu[2];
5244 m_fptotk = sqrt(1+pow(m_fhelixk[4],2))/m_fhelixk[2];
5245 m_fptotp = sqrt(1+pow(m_fhelixp[4],2))/m_fhelixp[2];
5247 m_zpt = 1/m_zhelix[2];
5248 m_zpte = 1/m_zhelixe[2];
5249 m_zptmu = 1/m_zhelixmu[2];
5250 m_zptk = 1/m_zhelixk[2];
5251 m_zptp = 1/m_zhelixp[2];
5253 m_fpt = 1/m_fhelix[2];
5254 m_fpte = 1/m_fhelixe[2];
5255 m_fptmu = 1/m_fhelixmu[2];
5256 m_fptk = 1/m_fhelixk[2];
5257 m_fptp = 1/m_fhelixp[2];
5259 m_lpt = 1/m_lhelix[2];
5260 m_lpte = 1/m_lhelixe[2];
5261 m_lptmu = 1/m_lhelixmu[2];
5262 m_lptk = 1/m_lhelixk[2];
5263 m_lptp = 1/m_lhelixp[2];
5265 m_lptot = sqrt(1+pow(m_lhelix[4],2))/m_lhelix[2];
5266 m_lptote = sqrt(1+pow(m_lhelixe[4],2))/m_lhelixe[2];
5267 m_lptotmu = sqrt(1+pow(m_lhelixmu[4],2))/m_lhelixmu[2];
5268 m_lptotk = sqrt(1+pow(m_lhelixk[4],2))/m_lhelixk[2];
5269 m_lptotp = sqrt(1+pow(m_lhelixp[4],2))/m_lhelixp[2];
5271 m_zptot = sqrt(1+pow(m_zhelix[4],2))/m_zhelix[2];
5272 m_zptote = sqrt(1+pow(m_zhelixe[4],2))/m_zhelixe[2];
5273 m_zptotmu = sqrt(1+pow(m_zhelixmu[4],2))/m_zhelixmu[2];
5274 m_zptotk = sqrt(1+pow(m_zhelixk[4],2))/m_zhelixk[2];
5275 m_zptotp = sqrt(1+pow(m_zhelixp[4],2))/m_zhelixp[2];
5277 m_zsigp = sqrt(pow((m_zptot/m_zhelix[2]),2)*m_zerror[5]+
5278 pow((m_zhelix[4]/m_zptot),2)*pow((1/m_zhelix[2]),4)*m_zerror[14]-
5279 2*m_zhelix[4]*m_zerror[12]*pow((1/m_zhelix[2]),3));
5280 m_zsigpe = sqrt(pow((m_zptote/m_zhelixe[2]),2)*m_zerrore[5]+
5281 pow((m_zhelixe[4]/m_zptote),2)*pow((1/m_zhelixe[2]),4)*m_zerrore[14]-
5282 2*m_zhelixe[4]*m_zerrore[12]*pow((1/m_zhelixe[2]),3));
5283 m_zsigpmu = sqrt(pow((m_zptotmu/m_zhelixmu[2]),2)*m_zerrormu[5]+
5284 pow((m_zhelixmu[4]/m_zptotmu),2)*pow((1/m_zhelixmu[2]),4)*m_zerrormu[14]-
5285 2*m_zhelixmu[4]*m_zerrormu[12]*pow((1/m_zhelixmu[2]),3));
5286 m_zsigpk = sqrt(pow((m_zptotk/m_zhelixk[2]),2)*m_zerrork[5]+
5287 pow((m_zhelixk[4]/m_zptotk),2)*pow((1/m_zhelixk[2]),4)*m_zerrork[14]-
5288 2*m_zhelixk[4]*m_zerrork[12]*pow((1/m_zhelixk[2]),3));
5289 m_zsigpp = sqrt(pow((m_zptotp/m_zhelixp[2]),2)*m_zerrorp[5]+
5290 pow((m_zhelixp[4]/m_zptotp),2)*pow((1/m_zhelixp[2]),4)*m_zerrorp[14]-
5291 2*m_zhelixp[4]*m_zerrorp[12]*pow((1/m_zhelixp[2]),3));
5294 StatusCode sc1 = m_nt1->write();
5295 if( sc1.isFailure() ) cout<<
"Ntuple1 filling failed!"<<endl;
5300 phi1 = (*iter_trk)->getFFi0();
5301 r1 = (*iter_trk)->getFDr();
5302 z1 = (*iter_trk)->getFDz();
5303 kap1 = (*iter_trk)->getFCpa();
5304 tanl1 = (*iter_trk)->getFTanl();
5307 p1 = sqrt(1+tanl1*tanl1)/kap1;
5308 the1 =
M_PI/2-atan(tanl1);
5309 }
else if(jj == 2) {
5310 phi2 = (*iter_trk)->getFFi0();
5311 r2 = (*iter_trk)->getFDr();
5312 z2 = (*iter_trk)->getFDz();
5313 kap2 = (*iter_trk)->getFCpa();
5314 tanl2 = (*iter_trk)->getFTanl();
5317 p2 = sqrt(1+tanl2*tanl2)/kap1;
5318 the2 =
M_PI/2-atan(tanl2);
5326 m_delthe = the1 + the2;
5329 StatusCode sc2 = m_nt2->write();
5330 if( sc2.isFailure() ) cout<<
"Ntuple2 filling failed!"<<endl;
5338 cout <<
"Kalfitting finished " << std::endl;