366 NTuplePtr nt1(
ntupleSvc(),
"FILE104/n101" );
368 if ( nt1 ) m_nt1 = nt1;
371 m_nt1 =
ntupleSvc()->book(
"FILE104/n101", CLID_ColumnWiseTuple,
"KalFit" );
375 status = m_nt1->addItem(
"trackid", m_trackid );
376 status = m_nt1->addItem(
"stat", 5, 2, m_stat );
377 status = m_nt1->addItem(
"ndf", 5, 2, m_ndf );
378 status = m_nt1->addItem(
"chisq", 5, 2, m_chisq );
379 status = m_nt1->addItem(
"length", 5, m_length );
380 status = m_nt1->addItem(
"tof", 5, m_tof );
381 status = m_nt1->addItem(
"nhits", 5, m_nhits );
382 status = m_nt1->addItem(
"zhelix", 5, m_zhelix );
383 status = m_nt1->addItem(
"zhelixe", 5, m_zhelixe );
384 status = m_nt1->addItem(
"zhelixmu", 5, m_zhelixmu );
385 status = m_nt1->addItem(
"zhelixk", 5, m_zhelixk );
386 status = m_nt1->addItem(
"zhelixp", 5, m_zhelixp );
387 status = m_nt1->addItem(
"zptot", m_zptot );
388 status = m_nt1->addItem(
"zptote", m_zptote );
389 status = m_nt1->addItem(
"zptotmu", m_zptotmu );
390 status = m_nt1->addItem(
"zptotk", m_zptotk );
391 status = m_nt1->addItem(
"zptotp", m_zptotp );
393 status = m_nt1->addItem(
"zpt", m_zpt );
394 status = m_nt1->addItem(
"zpte", m_zpte );
395 status = m_nt1->addItem(
"zptmu", m_zptmu );
396 status = m_nt1->addItem(
"zptk", m_zptk );
397 status = m_nt1->addItem(
"zptp", m_zptp );
399 status = m_nt1->addItem(
"fptot", m_fptot );
400 status = m_nt1->addItem(
"fptote", m_fptote );
401 status = m_nt1->addItem(
"fptotmu", m_fptotmu );
402 status = m_nt1->addItem(
"fptotk", m_fptotk );
403 status = m_nt1->addItem(
"fptotp", m_fptotp );
404 status = m_nt1->addItem(
"fpt", m_fpt );
405 status = m_nt1->addItem(
"fpte", m_fpte );
406 status = m_nt1->addItem(
"fptmu", m_fptmu );
407 status = m_nt1->addItem(
"fptk", m_fptk );
408 status = m_nt1->addItem(
"fptp", m_fptp );
410 status = m_nt1->addItem(
"lptot", m_lptot );
411 status = m_nt1->addItem(
"lptote", m_lptote );
412 status = m_nt1->addItem(
"lptotmu", m_lptotmu );
413 status = m_nt1->addItem(
"lptotk", m_lptotk );
414 status = m_nt1->addItem(
"lptotp", m_lptotp );
415 status = m_nt1->addItem(
"lpt", m_lpt );
416 status = m_nt1->addItem(
"lpte", m_lpte );
417 status = m_nt1->addItem(
"lptmu", m_lptmu );
418 status = m_nt1->addItem(
"lptk", m_lptk );
419 status = m_nt1->addItem(
"lptp", m_lptp );
421 status = m_nt1->addItem(
"zsigp", m_zsigp );
422 status = m_nt1->addItem(
"zsigpe", m_zsigpe );
423 status = m_nt1->addItem(
"zsigpmu", m_zsigpmu );
424 status = m_nt1->addItem(
"zsigpk", m_zsigpk );
425 status = m_nt1->addItem(
"zsigpp", m_zsigpp );
426 status = m_nt1->addItem(
"fhelix", 5, m_fhelix );
427 status = m_nt1->addItem(
"fhelixe", 5, m_fhelixe );
428 status = m_nt1->addItem(
"fhelixmu", 5, m_fhelixmu );
429 status = m_nt1->addItem(
"fhelixk", 5, m_fhelixk );
430 status = m_nt1->addItem(
"fhelixp", 5, m_fhelixp );
431 status = m_nt1->addItem(
"lhelix", 5, m_lhelix );
432 status = m_nt1->addItem(
"lhelixe", 5, m_lhelixe );
433 status = m_nt1->addItem(
"lhelixmu", 5, m_lhelixmu );
434 status = m_nt1->addItem(
"lhelixk", 5, m_lhelixk );
435 status = m_nt1->addItem(
"lhelixp", 5, m_lhelixp );
438 status = m_nt1->addItem(
"zerror", 15, m_zerror );
439 status = m_nt1->addItem(
"zerrore", 15, m_zerrore );
440 status = m_nt1->addItem(
"zerrormu", 15, m_zerrormu );
441 status = m_nt1->addItem(
"zerrork", 15, m_zerrork );
442 status = m_nt1->addItem(
"zerrorp", 15, m_zerrorp );
443 status = m_nt1->addItem(
"ferror", 15, m_ferror );
444 status = m_nt1->addItem(
"ferrore", 15, m_ferrore );
445 status = m_nt1->addItem(
"ferrormu", 15, m_ferrormu );
446 status = m_nt1->addItem(
"ferrork", 15, m_ferrork );
447 status = m_nt1->addItem(
"ferrorp", 15, m_ferrorp );
448 status = m_nt1->addItem(
"lerror", 15, m_lerror );
449 status = m_nt1->addItem(
"lerrore", 15, m_lerrore );
450 status = m_nt1->addItem(
"lerrormu", 15, m_lerrormu );
451 status = m_nt1->addItem(
"lerrork", 15, m_lerrork );
452 status = m_nt1->addItem(
"lerrorp", 15, m_lerrorp );
456 status = m_nt1->addItem(
"evtid", m_evtid );
457 status = m_nt1->addItem(
"mchelix", 5, m_mchelix );
458 status = m_nt1->addItem(
"mcptot", m_mcptot );
459 status = m_nt1->addItem(
"mcpid", m_mcpid );
461 if ( status.isFailure() ) cout <<
"Ntuple1 add item failed!" << endl;
468 NTuplePtr nt2(
ntupleSvc(),
"FILE104/n102" );
470 if ( nt2 ) m_nt2 = nt2;
473 m_nt2 =
ntupleSvc()->book(
"FILE104/n102", CLID_ColumnWiseTuple,
"KalFitComp" );
476 status2 = m_nt2->addItem(
"delx", m_delx );
477 status2 = m_nt2->addItem(
"dely", m_dely );
478 status2 = m_nt2->addItem(
"delz", m_delz );
479 status2 = m_nt2->addItem(
"delthe", m_delthe );
480 status2 = m_nt2->addItem(
"delphi", m_delphi );
481 status2 = m_nt2->addItem(
"delp", m_delp );
482 status2 = m_nt2->addItem(
"delpx", m_delpx );
483 status2 = m_nt2->addItem(
"delpy", m_delpy );
484 status2 = m_nt2->addItem(
"delpz", m_delpz );
486 if ( status2.isFailure() ) cout <<
"Ntuple2 add item failed!" << endl;
493 NTuplePtr nt3(
ntupleSvc(),
"FILE104/n103" );
495 if ( nt3 ) m_nt3 = nt3;
498 m_nt3 =
ntupleSvc()->book(
"FILE104/n103", CLID_ColumnWiseTuple,
"PatRec" );
501 status3 = m_nt3->addItem(
"trkhelix", 5, m_trkhelix );
502 status3 = m_nt3->addItem(
"trkptot", m_trkptot );
505 status3 = m_nt3->addItem(
"trkerror", 15, m_trkerror );
506 status3 = m_nt3->addItem(
"trksigp", m_trksigp );
508 status3 = m_nt3->addItem(
"trkndf", m_trkndf );
509 status3 = m_nt3->addItem(
"trkchisq", m_trkchisq );
510 if ( status3.isFailure() ) cout <<
"Ntuple3 add item failed!" << endl;
517 NTuplePtr nt4(
ntupleSvc(),
"FILE104/n104" );
519 if ( nt4 ) m_nt4 = nt4;
522 m_nt4 =
ntupleSvc()->book(
"FILE104/n104", CLID_ColumnWiseTuple,
"PatRecComp" );
525 status4 = m_nt4->addItem(
"trkdelx", m_trkdelx );
526 status4 = m_nt4->addItem(
"trkdely", m_trkdely );
527 status4 = m_nt4->addItem(
"trkdelz", m_trkdelz );
528 status4 = m_nt4->addItem(
"trkdelthe", m_trkdelthe );
529 status4 = m_nt4->addItem(
"trkdelphi", m_trkdelphi );
530 status4 = m_nt4->addItem(
"trkdelp", m_trkdelp );
531 if ( status4.isFailure() ) cout <<
"Ntuple4 add item failed!" << endl;
537 NTuplePtr nt5(
ntupleSvc(),
"FILE104/n105" );
539 if ( nt5 ) m_nt5 = nt5;
542 m_nt5 =
ntupleSvc()->book(
"FILE104/n105", CLID_ColumnWiseTuple,
"KalFitdChisq" );
545 status5 = m_nt5->addItem(
"dchi2", m_dchi2 );
546 status5 = m_nt5->addItem(
"masshyp", m_masshyp );
547 status5 = m_nt5->addItem(
"residual_estim", m_residest );
548 status5 = m_nt5->addItem(
"residual", m_residnew );
549 status5 = m_nt5->addItem(
"layer", m_layer );
550 status5 = m_nt5->addItem(
"kaldr", m_anal_dr );
551 status5 = m_nt5->addItem(
"kalphi0", m_anal_phi0 );
552 status5 = m_nt5->addItem(
"kalkappa", m_anal_kappa );
553 status5 = m_nt5->addItem(
"kaldz", m_anal_dz );
554 status5 = m_nt5->addItem(
"kaltanl", m_anal_tanl );
555 status5 = m_nt5->addItem(
"dr_ea", m_anal_ea_dr );
556 status5 = m_nt5->addItem(
"phi0_ea", m_anal_ea_phi0 );
557 status5 = m_nt5->addItem(
"kappa_ea", m_anal_ea_kappa );
558 status5 = m_nt5->addItem(
"dz_ea", m_anal_ea_dz );
559 status5 = m_nt5->addItem(
"tanl_ea", m_anal_ea_tanl );
560 if ( status5.isFailure() ) cout <<
"Ntuple5 add item failed!" << endl;
563 NTuplePtr nt6(
ntupleSvc(),
"FILE104/n106" );
565 if ( nt6 ) m_nt6 = nt6;
568 m_nt6 =
ntupleSvc()->book(
"FILE104/n106", CLID_ColumnWiseTuple,
"kal seg" );
571 status6 = m_nt6->addItem(
"docaInc", m_docaInc );
572 status6 = m_nt6->addItem(
"docaExc", m_docaExc );
573 status6 = m_nt6->addItem(
"tdr", m_tdrift );
574 status6 = m_nt6->addItem(
"layerid", m_layerid );
575 status6 = m_nt6->addItem(
"event", m_eventNo );
576 status6 = m_nt6->addItem(
"residualInc", m_residualInc );
577 status6 = m_nt6->addItem(
"residualExc", m_residualExc );
578 status6 = m_nt6->addItem(
"lr", m_lr );
579 status6 = m_nt6->addItem(
"dd", m_dd );
580 status6 = m_nt6->addItem(
"y", m_yposition );
582 if ( status6.isFailure() ) cout <<
"Ntuple6 add item failed!" << endl;
671 if ( sc.isFailure() )
673 error() <<
"beginRun failed" << endmsg;
674 return StatusCode::FAILURE;
679 MsgStream log(
msgSvc(), name() );
680 log << MSG::INFO <<
"in execute()" << endmsg;
722 StatusCode sc = service(
"MagneticFieldSvc", IMFSvc );
723 if ( sc != StatusCode::SUCCESS )
724 { log << MSG::ERROR <<
"Unable to open Magnetic field service" << endmsg; }
734 std::cout <<
" execute, referField from MagneticFieldSvc: "
742 IDataProviderSvc* evtSvc = NULL;
743 Gaudi::svcLocator()->service(
"EventDataSvc", evtSvc );
744 if ( evtSvc ) { log << MSG::INFO <<
"makeTds:event Svc has been found" << endmsg; }
747 log << MSG::FATAL <<
"makeTds:Could not find eventSvc" << endmsg;
748 return StatusCode::SUCCESS;
752 IDataManagerSvc* dataManSvc;
753 dataManSvc =
dynamic_cast<IDataManagerSvc*
>( evtSvc );
754 DataObject* aKalTrackCol;
755 evtSvc->findObject(
"/Event/Recon/RecMdcKalTrackCol", aKalTrackCol );
758 dataManSvc->clearSubTree(
"/Event/Recon/RecMdcKalTrackCol" );
759 evtSvc->unregisterObject(
"/Event/Recon/RecMdcKalTrackCol" );
762 kalsc = evtSvc->registerObject(
"/Event/Recon/RecMdcKalTrackCol", kalcol );
763 if ( kalsc.isFailure() )
765 log << MSG::FATAL <<
"Could not register RecMdcKalTrack" << endmsg;
766 return StatusCode::SUCCESS;
768 log << MSG::INFO <<
"RecMdcKalTrackCol registered successfully!" << endmsg;
770 DataObject* aKalHelixSegCol;
771 evtSvc->findObject(
"/Event/Recon/RecMdcKalHelixSegCol", aKalHelixSegCol );
772 if ( aKalHelixSegCol )
774 dataManSvc->clearSubTree(
"/Event/Recon/RecMdcKalHelixSegCol" );
775 evtSvc->unregisterObject(
"/Event/Recon/RecMdcKalHelixSegCol" );
778 kalsc = evtSvc->registerObject(
"/Event/Recon/RecMdcKalHelixSegCol", helixsegcol );
779 if ( kalsc.isFailure() )
781 log << MSG::FATAL <<
"Could not register RecMdcKalHelixSeg" << endmsg;
782 return StatusCode::SUCCESS;
784 log << MSG::INFO <<
"RecMdcKalHelixSegCol register successfully!" << endmsg;
795 { std::cout <<
"ERROR OCCUR when dynamic_cast in KalFitAlg::execute ...!!" << std::endl; }
797 SmartDataPtr<Event::EventHeader> eventHeader( eventSvc(),
"/Event/EventHeader" );
800 log << MSG::WARNING <<
"Could not find Event Header" << endmsg;
801 return StatusCode::FAILURE;
803 int eventNo = eventHeader->eventNumber();
804 int runNo = eventHeader->runNumber();
809 SmartDataPtr<RecEsTimeCol> estimeCol( eventSvc(),
"/Event/Recon/RecEsTimeCol" );
810 if ( estimeCol && estimeCol->size() )
812 RecEsTimeCol::iterator iter_evt = estimeCol->begin();
813 t0 = ( *iter_evt )->getTest();
818 log << MSG::WARNING <<
"Could not find EvTimeCol" << endmsg;
819 return StatusCode::SUCCESS;
823 { std::cout <<
"in KalFitAlg , we get the event start time = " << t0 << std::endl; }
826 SmartDataPtr<MdcDigiCol> mdcDigiCol( evtSvc,
"/Event/Digi/MdcDigiCol" );
827 if ( sc != StatusCode::SUCCESS )
829 log << MSG::FATAL <<
"Could not find MdcDigiCol!" << endmsg;
830 return StatusCode::SUCCESS;
842 m_evtid = eventHeader->eventNumber();
845 SmartDataPtr<McParticleCol> mcPartCol( eventSvc(),
"/Event/MC/McParticleCol" );
848 log << MSG::WARNING <<
"Could not find McParticle" << endmsg;
854 McParticleCol::iterator i_mcTrk = mcPartCol->begin();
855 for ( ; i_mcTrk != mcPartCol->end(); i_mcTrk++ )
857 if ( !( *i_mcTrk )->primaryParticle() )
continue;
858 const HepLorentzVector& mom( ( *i_mcTrk )->initialFourMomentum() );
859 const HepLorentzVector& pos = ( *i_mcTrk )->initialPosition();
860 log << MSG::DEBUG <<
"MCINFO:particleId=" << ( *i_mcTrk )->particleProperty()
861 <<
" theta=" << mom.theta() <<
" phi=" << mom.phi() <<
" px=" << mom.px()
862 <<
" py=" << mom.py() <<
" pz=" << mom.pz() << endmsg;
864 int pid = ( *i_mcTrk )->particleProperty();
865 if ( pid > 0 ) { charge = m_particleTable->particle( pid )->charge(); }
868 charge = m_particleTable->particle( -pid )->charge();
871 else { log << MSG::WARNING <<
"wrong particle id, please check data" << endmsg; }
873 Hep3Vector mom2( mom.px(), mom.py(), mom.pz() );
875 Helix mchelix( pos2, mom2, charge );
876 log << MSG::DEBUG <<
"charge of the track " << charge << endmsg;
877 if (
debug_ == 4 ) cout <<
"helix: " << mchelix.
a() << endl;
879 for (
int j = 0; j < 5; j++ ) { m_mchelix[j] = mchelix.
a()[j]; }
881 m_mcptot = sqrt( 1 + pow( m_mchelix[4], 2 ) ) / m_mchelix[2];
889 SmartDataPtr<RecMdcTrackCol> newtrkCol( eventSvc(),
"/Event/Recon/RecMdcTrackCol" );
892 log << MSG::FATAL <<
"Could not find RecMdcTrackCol" << endmsg;
893 return ( StatusCode::SUCCESS );
895 log << MSG::INFO <<
"Begin to make MdcRecTrkCol and MdcRecWirhitCol" << endmsg;
900 mtrkadd_mgr->clear();
904 double trkx1 = 0., trkx2 = 0., trky1 = 0., trky2 = 0., trkz1 = 0., trkz2 = 0., trkthe1 = 0., trkthe2 = 0., trkphi1 = 0., trkphi2 = 0., trkp1 = 0.,
905 trkp2 = 0., trkr1 = 0., trkr2 = 0., trkkap1 = 0., trkkap2 = 0., trktanl1 = 0., trktanl2 = 0.;
909 RecMdcTrackCol::iterator iter_trk = newtrkCol->begin();
910 for (
int kj = 1; iter_trk != newtrkCol->end(); iter_trk++, kj++ )
914 csmp3[kj - 1] = ( *iter_trk )->p3();
915 csmphi[kj - 1] = ( *iter_trk )->phi();
921 for (
int j = 0, ij = 0; j < 5; j++ )
923 m_trkhelix[j] = ( *iter_trk )->helix()[j];
926 for (
int k = 0; k <= j; k++, ij++ ) { m_trkerror[ij] = ( *iter_trk )->err()[j][k]; }
930 m_trkptot = sqrt( 1 + pow( m_trkhelix[4], 2 ) ) / m_trkhelix[2];
935 sqrt( pow( ( m_trkptot / m_trkhelix[2] ), 2 ) * m_trkerror[5] +
936 pow( ( m_trkhelix[4] / m_trkptot ), 2 ) * pow( ( 1 / m_trkhelix[2] ), 4 ) *
938 2 * m_trkhelix[4] * m_trkerror[12] * pow( ( 1 / m_trkhelix[2] ), 3 ) );
940 m_trkndf = ( *iter_trk )->ndof();
941 m_trkchisq = ( *iter_trk )->chi2();
943 if (
debug_ == 4 ) cout <<
"Ea from RecMdcTrackCol..." << ( *iter_trk )->err() << endl;
945 StatusCode sc3 = m_nt3->write();
946 if ( sc3.isFailure() ) cout <<
"Ntuple3 filling failed!" << endl;
978 log << MSG::DEBUG <<
"retrieved MDC tracks:"
979 <<
" Nhits " << ( *iter_trk )->getNhits() <<
" Nster " << ( *iter_trk )->nster()
982 HitRefVec gothits = ( *iter_trk )->getVecHits();
986 rectrk->
id = ( *iter_trk )->trackId();
987 rectrk->
chiSq = ( *iter_trk )->chi2();
988 rectrk->
ndf = ( *iter_trk )->ndof();
989 rectrk->
fiTerm = ( *iter_trk )->getFiTerm();
990 rectrk->
nhits = ( *iter_trk )->getNhits();
991 rectrk->
nster = ( *iter_trk )->nster();
993 rectrk->
stat = ( *iter_trk )->stat();
994 status_temp = ( *iter_trk )->stat();
996 trkadd->
id = ( *iter_trk )->trackId();
1000 trkadd->
body = rectrk;
1001 rectrk->
add = trkadd;
1003 for (
int i = 0; i < 5; i++ )
1005 rectrk->
helix[i] = ( *iter_trk )->helix()[i];
1006 if ( i < 3 ) rectrk->
pivot[i] = ( *iter_trk )->getPivot()[i];
1007 for (
int j = 1; j < i + 2; j++ )
1008 { rectrk->
error[i * ( i + 1 ) / 2 + j - 1] = ( *iter_trk )->err()( i + 1, j ); }
1011 std::sort( gothits.begin(), gothits.end(), order_rechits );
1013 HitRefVec::iterator it_gothit = gothits.begin();
1014 for ( ; it_gothit != gothits.end(); it_gothit++ )
1017 if ( ( *it_gothit )->getStat() != 1 )
1021 log << MSG::WARNING <<
"this hit is not used in helix fitting!" << endmsg;
1026 log << MSG::DEBUG <<
"retrieved hits in MDC tracks:"
1027 <<
" hits DDL " << ( *it_gothit )->getDriftDistLeft() <<
" hits DDR "
1028 << ( *it_gothit )->getDriftDistRight() <<
" error DDL "
1029 << ( *it_gothit )->getErrDriftDistLeft() <<
" error DDR "
1030 << ( *it_gothit )->getErrDriftDistRight() <<
" id of hit "
1031 << ( *it_gothit )->getId() <<
" track id of hit " << ( *it_gothit )->getTrkId()
1032 <<
" hits ADC " << ( *it_gothit )->getAdc() << endmsg;
1035 whit->
id = ( *it_gothit )->getId();
1036 whit->
ddl = ( *it_gothit )->getDriftDistLeft();
1037 whit->
ddr = ( *it_gothit )->getDriftDistRight();
1038 whit->
erddl = ( *it_gothit )->getErrDriftDistLeft();
1039 whit->
erddr = ( *it_gothit )->getErrDriftDistRight();
1040 whit->
pChiSq = ( *it_gothit )->getChisqAdd();
1041 whit->
lr = ( *it_gothit )->getFlagLR();
1042 whit->
stat = ( *it_gothit )->getStat();
1043 mdcid = ( *it_gothit )->getMdcId();
1047 int wid = w0id + localwid;
1048 log << MSG::INFO <<
"lr from PR: " << whit->
lr <<
" layerId = " << layid
1049 <<
" wireId = " << localwid << endmsg;
1058 whit->
tdc = ( *it_gothit )->getTdc();
1059 whit->
adc = ( *it_gothit )->getAdc();
1060 rectrk->
hitcol.push_back( whit );
1061 mhit_mgr->push_back( *whit );
1063 mtrk_mgr->push_back( *rectrk );
1064 mtrkadd_mgr->push_back( *trkadd );
1073 m_trkdelx = trkx1 - trkx2;
1074 m_trkdely = trky1 - trky2;
1075 m_trkdelz = trkz1 - trkz2;
1076 m_trkdelthe = trkthe1 + trkthe2;
1077 m_trkdelphi = trkphi1 - trkphi2;
1078 m_trkdelp = trkp1 - trkp2;
1079 StatusCode sc4 = m_nt4->write();
1080 if ( sc4.isFailure() ) cout <<
"Ntuple4 filling failed!" << endl;
1085 std::cout <<
"before refit,ntrk,nhits,nadd..." << mtrk_mgr->size() <<
"********"
1086 << mhit_mgr->size() <<
"****" << mtrkadd_mgr->size() << endl;
1092 double mdang = 180.0 - csmp3[0].angle( csmp3[1].
unit() ) * 180.0 /
M_PI;
1093 double mdphi = 180.0 - fabs( csmphi[0] - csmphi[1] ) * 180.0 /
M_PI;
1097 if (
usage_ == 2 && ( mtrk_mgr->size() ) == 2 && fabs( mdang ) <
m_dangcut &&
1100 if (
usage_ == 3 && ( mtrk_mgr->size() ) == 1 && status_temp == -1 )
1103 log << MSG::DEBUG <<
"after kalman_fitting(),but in execute...." << endmsg;
1109 SmartDataPtr<RecMdcKalTrackCol> recmdckaltrkCol( eventSvc(),
1110 "/Event/Recon/RecMdcKalTrackCol" );
1115 RecMdcKalTrackCol::iterator KalTrk = recmdckaltrkCol->begin();
1117 for ( ; KalTrk != recmdckaltrkCol->end(); KalTrk++ )
1120 for (
int hypo = 0; hypo < 5; hypo++ )
1122 if ( ( *KalTrk )->getStat( 0, hypo ) == 1 ) nFailedTrks[hypo]++;
1125 HelixSegRefVec::iterator iter_hit = gothelixsegs.begin();
1127 int nhitofthistrk = 0;
1128 for ( ; iter_hit != gothelixsegs.end(); iter_hit++ )
1136 iter_hit = gothelixsegs.begin();
1178 return StatusCode::SUCCESS;
3120 MsgStream log(
msgSvc(), name() );
3121 double Pt_threshold( 0.3 );
3122 Hep3Vector IP( 0, 0, 0 );
3128 if ( !&whMgr )
return;
3131 int ntrk = mdcMgr->size();
3132 double* rPt =
new double[ntrk];
3133 int* rOM =
new int[ntrk];
3134 unsigned int* order =
new unsigned int[ntrk];
3135 unsigned int* rCont =
new unsigned int[ntrk];
3136 unsigned int* rGen =
new unsigned int[ntrk];
3139 for ( vector<MdcRec_trk>::iterator it = mdcMgr->begin(), end = mdcMgr->end(); it != end;
3144 if ( it->helix[2] ) rPt[index] = 1 / fabs( it->helix[2] );
3145 if (
debug_ == 4 ) cout <<
"rPt...[ " << index <<
" ]...." << rPt[index] << endl;
3146 if ( rPt[index] < 0 ) rPt[index] =
DBL_MAX;
3148 std::vector<MdcRec_wirhit*> pt = it->hitcol;
3149 if (
debug_ == 4 ) cout <<
"ppt size: " << pt.size() << endl;
3150 int outermost( -1 );
3151 for ( vector<MdcRec_wirhit*>::iterator ii = pt.end() - 1; ii != pt.begin() - 1; ii-- )
3153 int lyr( ( *ii )->geo->Lyr()->Id() );
3154 if ( outermost < lyr ) outermost = lyr;
3155 if (
debug_ == 4 ) cout <<
"outmost: " << outermost <<
" lyr: " << lyr << endl;
3157 rOM[index] = outermost;
3158 order[index] = index;
3163 for (
int j, k = ntrk - 1; k >= 0; k = j )
3166 for (
int i = 1; i <= k; i++ )
3167 if ( rPt[i - 1] < rPt[i] )
3170 std::swap( order[i], order[j] );
3171 std::swap( rPt[i], rPt[j] );
3172 std::swap( rOM[i], rOM[j] );
3173 std::swap( rCont[i], rCont[j] );
3174 std::swap( rGen[i], rGen[j] );
3182 DataObject* aReconEvent;
3183 eventSvc()->findObject(
"/Event/Recon", aReconEvent );
3188 StatusCode sc = eventSvc()->registerObject(
"/Event/Recon", recevt );
3189 if ( sc != StatusCode::SUCCESS )
3191 log << MSG::FATAL <<
"Could not register ReconEvent" << endmsg;
3199 log << MSG::INFO <<
"beginning to make RecMdcKalTrackCol" << endmsg;
3201 for (
int l = 0; l < ntrk; l++ )
3212 MdcRec_trk& TrasanTRK = *( mdcMgr->begin() + order[l] );
3213 MdcRec_trk_add& TrasanTRK_add = *( mdc_addMgr->begin() + order[l] );
3216 int trasqual = TrasanTRK_add.
quality;
3217 if (
debug_ == 4 ) cout <<
"kalman_fitting>trasqual... " << trasqual << endl;
3218 if ( trasqual )
continue;
3221 if (
debug_ == 4 ) cout <<
"******* KalFit NUMBER : " << newcount << std::endl;
3225 if ( ( TrasanTRK_add.
decision & 32 ) == 32 || ( TrasanTRK_add.
decision & 64 ) == 64 )
3232 for (
int i = 0; i < 5; i++ ) a[i] = TrasanTRK.
helix[i];
3234 HepSymMatrix ea( 5 );
3235 for (
int i = 0, k = 0; i < 5; i++ )
3237 for (
int j = 0; j <= i; j++ )
3240 ea[j][i] = ea[i][j];
3243 double fiTerm = TrasanTRK.
fiTerm;
3251 int inlyr( 999 ), outlyr( -1 );
3252 int* rStat =
new int[43];
3253 for (
int irStat = 0; irStat < 43; ++irStat ) rStat[irStat] = 0;
3254 std::vector<MdcRec_wirhit*> pt = TrasanTRK.
hitcol;
3256 if (
debug_ == 4 ) cout <<
"*********Pt size****" << pt.size() << endl;
3258 int Num[43] = { 0 };
3259 for ( vector<MdcRec_wirhit*>::iterator ii = pt.end() - 1; ii != pt.begin() - 1; ii-- )
3260 { Num[( *ii )->geo->Lyr()->Id()]++; }
3262 for ( vector<MdcRec_wirhit*>::iterator ii = pt.end() - 1; ii != pt.begin() - 1; ii-- )
3265 if ( Num[( *ii )->geo->Lyr()->Id()] > 3 )
3268 cout <<
"WARNING: I found " << Num[( *ii )->geo->Lyr()->Id()]
3269 <<
" hits in the layer " << ( *ii )->geo->Lyr()->Id() << std::endl;
3289 double dist[2] = { rechit.
ddl, rechit.
ddr };
3290 double erdist[2] = { rechit.
erddl, rechit.
erddr };
3293 int lr_decision( 0 );
3296 if ( rechit.
lr == 2 || rechit.
lr == 0 ) lr_decision = -1;
3298 else if ( rechit.
lr == 1 ) lr_decision = 1;
3301 int ind( geo->
Lyr()->
Id() );
3308 if ( inlyr > ind ) inlyr = ind;
3309 if ( outlyr < ind ) outlyr = ind;
3311 if (
debug_ == 4 ) cout <<
"**** NUMBER OF Mdc HITS (TRASAN) = " << hit_asso << std::endl;
3314 int empty_between( 0 ), empty( 0 );
3315 for (
int i = inlyr; i <= outlyr; i++ )
3316 if ( !rStat[i] ) empty_between++;
3317 empty = empty_between + inlyr + ( 42 - outlyr );
3324 for ( std::vector<KalFitHitMdc>::iterator it_hit = track_lead.
HitsMdc().begin();
3325 it_hit != track_lead.
HitsMdc().end(); it_hit++ )
3334 unsigned int nhit = track_lead.
HitsMdc().size();
3335 if ( !nhit &&
debug_ == 4 )
3337 cout <<
" ATTENTION TRACK WITH ONLY HITS " << nhit << std::endl;
3342 double KalFitst( 0 ), KalFitax( 0 ), KalFitschi2( 0 );
3345 Hep3Vector outer_pivot( track_lead.
x( fiTerm ) );
3347 track_lead.
pivot( outer_pivot );
3352 HepSymMatrix Eakal( 5, 0 );
3355 double costheta = track_lead.
a()[4] / sqrt( 1.0 + track_lead.
a()[4] * track_lead.
a()[4] );
3356 if ( ( 1.0 / fabs( track_lead.
a()[2] ) <
pt_cut_ ) && ( fabs( costheta ) >
theta_cut_ ) )
3364 cout <<
"from Mdc Pattern Recognition: " << std::endl;
3366 Helix work( track_lead.
pivot(), track_lead.
a(), track_lead.
Ea() );
3368 cout <<
" dr = " << work.
a()[0] <<
", Er_dr = " << sqrt( work.
Ea()[0][0] ) << std::endl;
3369 cout <<
" phi0 = " << work.
a()[1] <<
", Er_phi0 = " << sqrt( work.
Ea()[1][1] )
3371 cout <<
" PT = " << 1 / work.
a()[2] <<
", Er_kappa = " << sqrt( work.
Ea()[2][2] )
3373 cout <<
" dz = " << work.
a()[3] <<
", Er_dz = " << sqrt( work.
Ea()[3][3] ) << std::endl;
3374 cout <<
" tanl = " << work.
a()[4] <<
", Er_tanl = " << sqrt( work.
Ea()[4][4] )
3386 cout <<
" Mdc FIRST KALMAN FIT " << std::endl;
3387 Helix work( track_lead.
pivot(), track_lead.
a(), track_lead.
Ea() );
3389 cout <<
" dr = " << work.
a()[0] <<
", Er_dr = " << sqrt( work.
Ea()[0][0] ) << std::endl;
3390 cout <<
" phi0 = " << work.
a()[1] <<
", Er_phi0 = " << sqrt( work.
Ea()[1][1] )
3392 cout <<
" PT = " << 1 / work.
a()[2] <<
", Er_kappa = " << sqrt( work.
Ea()[2][2] )
3394 cout <<
" dz = " << work.
a()[3] <<
", Er_dz = " << sqrt( work.
Ea()[3][3] ) << std::endl;
3395 cout <<
" tanl = " << work.
a()[4] <<
", Er_tanl = " << sqrt( work.
Ea()[4][4] )
3402 complete_track( TrasanTRK, TrasanTRK_add, track_lead, kaltrk, kalcol, segcol, 1 );
3405 IDataProviderSvc* eventSvc = NULL;
3406 Gaudi::svcLocator()->service(
"EventDataSvc", eventSvc );
3407 if ( eventSvc ) { log << MSG::INFO <<
"makeTds:event Svc has been found" << endmsg; }
3410 log << MSG::FATAL <<
"makeTds:Could not find eventSvc" << endmsg;
3415 IDataManagerSvc* dataManSvc;
3416 dataManSvc =
dynamic_cast<IDataManagerSvc*
>( eventSvc );
3417 DataObject* aKalTrackCol;
3418 eventSvc->findObject(
"/Event/Recon/RecMdcKalTrackCol", aKalTrackCol );
3421 dataManSvc->clearSubTree(
"/Event/Recon/RecMdcKalTrackCol" );
3422 eventSvc->unregisterObject(
"/Event/Recon/RecMdcKalTrackCol" );
3425 kalsc = eventSvc->registerObject(
"/Event/Recon/RecMdcKalTrackCol", kalcol );
3426 if ( kalsc.isFailure() )
3428 log << MSG::FATAL <<
"Could not register RecMdcKalTrack" << endmsg;
3431 log << MSG::INFO <<
"RecMdcKalTrackCol registered successfully!" << endmsg;
3435 DataObject* aRecKalSegEvent;
3436 eventSvc->findObject(
"/Event/Recon/RecMdcKalHelixSegCol", aRecKalSegEvent );
3437 if ( aRecKalSegEvent )
3440 segsc = eventSvc->unregisterObject(
"/Event/Recon/RecMdcKalHelixSegCol" );
3441 if ( segsc != StatusCode::SUCCESS )
3443 log << MSG::FATAL <<
"Could not unregister RecMdcKalHelixSegCol collection" << endmsg;
3448 segsc = eventSvc->registerObject(
"/Event/Recon/RecMdcKalHelixSegCol", segcol );
3449 if ( segsc.isFailure() )
3451 log << MSG::FATAL <<
"Could not register RecMdcKalHelixSeg" << endmsg;
3454 log << MSG::INFO <<
"RecMdcKalHelixSegCol registered successfully!" << endmsg;
3456 double x1( 0. ), x2( 0. ), y1( 0. ), y2( 0. ), z1( 0. ), z2( 0. ), the1( 0. ), the2( 0. ),
3458 double r1( 0. ), r2( 0. ), kap1( 999. ), kap2( 999. ), tanl1( 0 ), tanl2( 0. );
3459 double px1( 0. ), px2( 0. ), py1( 0. ), py2( 0. ), pz1( 0. ), pz2( 0. ), charge1( 0. ),
3463 SmartDataPtr<RecMdcKalTrackCol> kaltrkCol( eventSvc,
"/Event/Recon/RecMdcKalTrackCol" );
3466 log << MSG::FATAL <<
"Could not find RecMdcKalTrackCol" << endmsg;
3469 log << MSG::INFO <<
"Begin to check RecMdcKalTrackCol" << endmsg;
3470 RecMdcKalTrackCol::iterator iter_trk = kaltrkCol->begin();
3471 for (
int jj = 1; iter_trk != kaltrkCol->end(); iter_trk++, jj++ )
3473 log << MSG::DEBUG <<
"retrieved MDC Kalmantrack:"
3474 <<
"Track Id: " << ( *iter_trk )->getTrackId()
3475 <<
" Mass of the fit: " << ( *iter_trk )->getMass( 2 ) << endmsg
3476 <<
" Length of the track: " << ( *iter_trk )->getLength( 2 )
3477 <<
" Tof of the track: " << ( *iter_trk )->getTof( 2 ) << endmsg
3478 <<
" Chisq of the fit: " << ( *iter_trk )->getChisq( 0, 2 ) <<
" "
3479 << ( *iter_trk )->getChisq( 1, 2 ) << endmsg
3480 <<
"Ndf of the fit: " << ( *iter_trk )->getNdf( 0, 2 ) <<
" "
3481 << ( *iter_trk )->getNdf( 1, 2 ) << endmsg <<
"Kappa " << ( *iter_trk )->getZHelix()[2]
3483 for (
int i = 0; i < 43; i++ )
3485 log << MSG::DEBUG <<
"retrieved pathl[" << i <<
"]= " << ( *iter_trk )->getPathl( i )
3491 m_trackid = ( *iter_trk )->getTrackId();
3493 for (
int jj = 0, iii = 0; jj < 5; jj++ )
3496 m_length[jj] = ( *iter_trk )->getLength( jj );
3497 m_tof[jj] = ( *iter_trk )->getTof( jj );
3498 m_nhits[jj] = ( *iter_trk )->getNhits( jj );
3499 m_zhelix[jj] = ( *iter_trk )->getZHelix()[jj];
3500 m_zhelixe[jj] = ( *iter_trk )->getZHelixE()[jj];
3501 m_zhelixmu[jj] = ( *iter_trk )->getZHelixMu()[jj];
3502 m_zhelixk[jj] = ( *iter_trk )->getZHelixK()[jj];
3503 m_zhelixp[jj] = ( *iter_trk )->getZHelixP()[jj];
3504 m_fhelix[jj] = ( *iter_trk )->getFHelix()[jj];
3505 m_fhelixe[jj] = ( *iter_trk )->getFHelixE()[jj];
3506 m_fhelixmu[jj] = ( *iter_trk )->getFHelixMu()[jj];
3507 m_fhelixk[jj] = ( *iter_trk )->getFHelixK()[jj];
3508 m_fhelixp[jj] = ( *iter_trk )->getFHelixP()[jj];
3509 m_lhelix[jj] = ( *iter_trk )->getLHelix()[jj];
3510 m_lhelixe[jj] = ( *iter_trk )->getLHelixE()[jj];
3511 m_lhelixmu[jj] = ( *iter_trk )->getLHelixMu()[jj];
3512 m_lhelixk[jj] = ( *iter_trk )->getLHelixK()[jj];
3513 m_lhelixp[jj] = ( *iter_trk )->getLHelixP()[jj];
3517 for (
int kk = 0; kk <= jj; kk++, iii++ )
3519 m_zerror[iii] = ( *iter_trk )->getZError()[jj][kk];
3520 m_zerrore[iii] = ( *iter_trk )->getZErrorE()[jj][kk];
3521 m_zerrormu[iii] = ( *iter_trk )->getZErrorMu()[jj][kk];
3522 m_zerrork[iii] = ( *iter_trk )->getZErrorK()[jj][kk];
3523 m_zerrorp[iii] = ( *iter_trk )->getZErrorP()[jj][kk];
3524 m_ferror[iii] = ( *iter_trk )->getFError()[jj][kk];
3525 m_ferrore[iii] = ( *iter_trk )->getFErrorE()[jj][kk];
3526 m_ferrormu[iii] = ( *iter_trk )->getFErrorMu()[jj][kk];
3527 m_ferrork[iii] = ( *iter_trk )->getFErrorK()[jj][kk];
3528 m_ferrorp[iii] = ( *iter_trk )->getFErrorP()[jj][kk];
3529 m_lerror[iii] = ( *iter_trk )->getLError()[jj][kk];
3530 m_lerrore[iii] = ( *iter_trk )->getLErrorE()[jj][kk];
3531 m_lerrormu[iii] = ( *iter_trk )->getLErrorMu()[jj][kk];
3532 m_lerrork[iii] = ( *iter_trk )->getLErrorK()[jj][kk];
3533 m_lerrorp[iii] = ( *iter_trk )->getLErrorP()[jj][kk];
3572 m_chisq[0][0] = ( *iter_trk )->getChisq( 0, 0 );
3573 m_chisq[1][0] = ( *iter_trk )->getChisq( 0, 1 );
3574 m_chisq[2][0] = ( *iter_trk )->getChisq( 0, 2 );
3575 m_chisq[3][0] = ( *iter_trk )->getChisq( 0, 3 );
3576 m_chisq[4][0] = ( *iter_trk )->getChisq( 0, 4 );
3577 m_chisq[0][1] = ( *iter_trk )->getChisq( 1, 0 );
3578 m_chisq[1][1] = ( *iter_trk )->getChisq( 1, 1 );
3579 m_chisq[2][1] = ( *iter_trk )->getChisq( 1, 2 );
3580 m_chisq[3][1] = ( *iter_trk )->getChisq( 1, 3 );
3581 m_chisq[4][1] = ( *iter_trk )->getChisq( 1, 4 );
3583 m_ndf[0][0] = ( *iter_trk )->getNdf( 0, 0 );
3584 m_ndf[1][0] = ( *iter_trk )->getNdf( 0, 1 );
3585 m_ndf[2][0] = ( *iter_trk )->getNdf( 0, 2 );
3586 m_ndf[3][0] = ( *iter_trk )->getNdf( 0, 3 );
3587 m_ndf[4][0] = ( *iter_trk )->getNdf( 0, 4 );
3588 m_ndf[0][1] = ( *iter_trk )->getNdf( 1, 0 );
3589 m_ndf[1][1] = ( *iter_trk )->getNdf( 1, 1 );
3590 m_ndf[2][1] = ( *iter_trk )->getNdf( 1, 2 );
3591 m_ndf[3][1] = ( *iter_trk )->getNdf( 1, 3 );
3592 m_ndf[4][1] = ( *iter_trk )->getNdf( 1, 4 );
3594 m_stat[0][0] = ( *iter_trk )->getStat( 0, 0 );
3595 m_stat[1][0] = ( *iter_trk )->getStat( 0, 1 );
3596 m_stat[2][0] = ( *iter_trk )->getStat( 0, 2 );
3597 m_stat[3][0] = ( *iter_trk )->getStat( 0, 3 );
3598 m_stat[4][0] = ( *iter_trk )->getStat( 0, 4 );
3599 m_stat[0][1] = ( *iter_trk )->getStat( 1, 0 );
3600 m_stat[1][1] = ( *iter_trk )->getStat( 1, 1 );
3601 m_stat[2][1] = ( *iter_trk )->getStat( 1, 2 );
3602 m_stat[3][1] = ( *iter_trk )->getStat( 1, 3 );
3603 m_stat[4][1] = ( *iter_trk )->getStat( 1, 4 );
3605 m_fptot = sqrt( 1 + pow( m_fhelix[4], 2 ) ) / m_fhelix[2];
3606 m_fptote = sqrt( 1 + pow( m_fhelixe[4], 2 ) ) / m_fhelixe[2];
3607 m_fptotmu = sqrt( 1 + pow( m_fhelixmu[4], 2 ) ) / m_fhelixmu[2];
3608 m_fptotk = sqrt( 1 + pow( m_fhelixk[4], 2 ) ) / m_fhelixk[2];
3609 m_fptotp = sqrt( 1 + pow( m_fhelixp[4], 2 ) ) / m_fhelixp[2];
3611 m_lptot = sqrt( 1 + pow( m_lhelix[4], 2 ) ) / m_lhelix[2];
3612 m_lptote = sqrt( 1 + pow( m_lhelixe[4], 2 ) ) / m_lhelixe[2];
3613 m_lptotmu = sqrt( 1 + pow( m_lhelixmu[4], 2 ) ) / m_lhelixmu[2];
3614 m_lptotk = sqrt( 1 + pow( m_lhelixk[4], 2 ) ) / m_lhelixk[2];
3615 m_lptotp = sqrt( 1 + pow( m_lhelixp[4], 2 ) ) / m_lhelixp[2];
3617 m_lpt = 1 / m_lhelix[2];
3618 m_lpte = 1 / m_lhelixe[2];
3619 m_lptmu = 1 / m_lhelixmu[2];
3620 m_lptk = 1 / m_lhelixk[2];
3621 m_lptp = 1 / m_lhelixp[2];
3623 m_fpt = 1 / m_fhelix[2];
3624 m_fpte = 1 / m_fhelixe[2];
3625 m_fptmu = 1 / m_fhelixmu[2];
3626 m_fptk = 1 / m_fhelixk[2];
3627 m_fptp = 1 / m_fhelixp[2];
3631 std::cout <<
" " << std::endl;
3632 std::cout <<
"in file Kalman_fitting_anal ,the m_fpt is .." << m_fpt << std::endl;
3633 std::cout <<
"in file Kalman_fitting_anal ,the m_fpte is .." << m_fpte << std::endl;
3634 std::cout <<
"in file Kalman_fitting_anal ,the m_fptmu is .." << m_fptmu << std::endl;
3635 std::cout <<
"in file Kalman_fitting_anal ,the m_fptk is .." << m_fptk << std::endl;
3636 std::cout <<
"in file Kalman_fitting_anal ,the m_fptp is .." << m_fptp << std::endl;
3639 m_zpt = 1 / m_zhelix[2];
3640 m_zpte = 1 / m_zhelixe[2];
3641 m_zptmu = 1 / m_zhelixmu[2];
3642 m_zptk = 1 / m_zhelixk[2];
3643 m_zptp = 1 / m_zhelixp[2];
3647 std::cout <<
"in file Kalman_fitting_anal ,the m_zpt is .." << m_zpt << std::endl;
3648 std::cout <<
"in file Kalman_fitting_anal ,the m_zpte is .." << m_zpte << std::endl;
3649 std::cout <<
"in file Kalman_fitting_anal ,the m_zptmu is .." << m_zptmu << std::endl;
3650 std::cout <<
"in file Kalman_fitting_anal ,the m_zptk is .." << m_zptk << std::endl;
3651 std::cout <<
"in file Kalman_fitting_anal ,the m_zptp is .." << m_zptp << std::endl;
3653 m_zptot = sqrt( 1 + pow( m_zhelix[4], 2 ) ) / m_zhelix[2];
3654 m_zptote = sqrt( 1 + pow( m_zhelixe[4], 2 ) ) / m_zhelixe[2];
3655 m_zptotmu = sqrt( 1 + pow( m_zhelixmu[4], 2 ) ) / m_zhelixmu[2];
3656 m_zptotk = sqrt( 1 + pow( m_zhelixk[4], 2 ) ) / m_zhelixk[2];
3657 m_zptotp = sqrt( 1 + pow( m_zhelixp[4], 2 ) ) / m_zhelixp[2];
3661 std::cout <<
"in file Kalman_fitting_anal ,the m_zptot is .." << m_zptot << std::endl;
3662 std::cout <<
"in file Kalman_fitting_anal ,the m_zptote is .." << m_zptote
3664 std::cout <<
"in file Kalman_fitting_anal ,the m_zptotmu is .." << m_zptotmu
3666 std::cout <<
"in file Kalman_fitting_anal ,the m_zptotk is .." << m_zptotk
3668 std::cout <<
"in file Kalman_fitting_anal ,the m_zptotp is .." << m_zptotp
3674 m_zsigp = sqrt( pow( ( m_zptot / m_zhelix[2] ), 2 ) * m_zerror[5] +
3675 pow( ( m_zhelix[4] / m_zptot ), 2 ) * pow( ( 1 / m_zhelix[2] ), 4 ) *
3677 2 * m_zhelix[4] * m_zerror[12] * pow( ( 1 / m_zhelix[2] ), 3 ) );
3678 m_zsigpe = sqrt( pow( ( m_zptote / m_zhelixe[2] ), 2 ) * m_zerrore[5] +
3679 pow( ( m_zhelixe[4] / m_zptote ), 2 ) *
3680 pow( ( 1 / m_zhelixe[2] ), 4 ) * m_zerrore[14] -
3681 2 * m_zhelixe[4] * m_zerrore[12] * pow( ( 1 / m_zhelixe[2] ), 3 ) );
3683 sqrt( pow( ( m_zptotmu / m_zhelixmu[2] ), 2 ) * m_zerrormu[5] +
3684 pow( ( m_zhelixmu[4] / m_zptotmu ), 2 ) * pow( ( 1 / m_zhelixmu[2] ), 4 ) *
3686 2 * m_zhelixmu[4] * m_zerrormu[12] * pow( ( 1 / m_zhelixmu[2] ), 3 ) );
3687 m_zsigpk = sqrt( pow( ( m_zptotk / m_zhelixk[2] ), 2 ) * m_zerrork[5] +
3688 pow( ( m_zhelixk[4] / m_zptotk ), 2 ) *
3689 pow( ( 1 / m_zhelixk[2] ), 4 ) * m_zerrork[14] -
3690 2 * m_zhelixk[4] * m_zerrork[12] * pow( ( 1 / m_zhelixk[2] ), 3 ) );
3691 m_zsigpp = sqrt( pow( ( m_zptotp / m_zhelixp[2] ), 2 ) * m_zerrorp[5] +
3692 pow( ( m_zhelixp[4] / m_zptotp ), 2 ) *
3693 pow( ( 1 / m_zhelixp[2] ), 4 ) * m_zerrorp[14] -
3694 2 * m_zhelixp[4] * m_zerrorp[12] * pow( ( 1 / m_zhelixp[2] ), 3 ) );
3697 StatusCode sc1 = m_nt1->write();
3698 if ( sc1.isFailure() ) cout <<
"Ntuple1 filling failed!" << endl;
3705 phi1 = ( *iter_trk )->getFFi0();
3706 r1 = ( *iter_trk )->getFDr();
3707 z1 = ( *iter_trk )->getFDz();
3708 kap1 = ( *iter_trk )->getFCpa();
3709 tanl1 = ( *iter_trk )->getFTanl();
3710 charge1 = kap1 / fabs( kap1 );
3713 p1 = sqrt( 1 + tanl1 * tanl1 ) / kap1;
3714 the1 =
M_PI / 2 - atan( tanl1 );
3715 px1 = -
sin(
phi1 ) / fabs( kap1 );
3716 py1 =
cos(
phi1 ) / fabs( kap1 );
3717 pz1 = tanl1 / fabs( kap1 );
3721 phi2 = ( *iter_trk )->getFFi0();
3722 r2 = ( *iter_trk )->getFDr();
3723 z2 = ( *iter_trk )->getFDz();
3724 kap2 = ( *iter_trk )->getFCpa();
3725 tanl2 = ( *iter_trk )->getFTanl();
3726 charge2 = kap2 / fabs( kap2 );
3729 p2 = sqrt( 1 + tanl2 * tanl2 ) / kap1;
3730 the2 =
M_PI / 2 - atan( tanl2 );
3731 px2 = -
sin(
phi2 ) / fabs( kap2 );
3732 py2 =
cos(
phi2 ) / fabs( kap2 );
3733 pz2 = tanl2 / fabs( kap2 );
3742 m_delthe = the1 + the2;
3745 m_delpx = charge1 * fabs( px1 ) + charge2 * fabs( px2 );
3746 m_delpy = charge1 * fabs( py1 ) + charge2 * fabs( py2 );
3747 m_delpz = charge1 * fabs( pz1 ) + charge2 * fabs( pz2 );
3749 StatusCode sc2 = m_nt2->write();
3750 if ( sc2.isFailure() ) cout <<
"Ntuple2 filling failed!" << endl;
3758 if (
debug_ == 4 ) cout <<
"Kalfitting finished " << std::endl;
3763 MsgStream log(
msgSvc(), name() );
3764 double Pt_threshold( 0.3 );
3765 Hep3Vector IP( 0, 0, 0 );
3772 if ( !&whMgr )
return;
3775 int ntrk = mdcMgr->size();
3776 double* rPt =
new double[ntrk];
3777 int* rOM =
new int[ntrk];
3778 unsigned int* order =
new unsigned int[ntrk];
3779 unsigned int* rCont =
new unsigned int[ntrk];
3780 unsigned int* rGen =
new unsigned int[ntrk];
3783 for ( vector<MdcRec_trk>::iterator it = mdcMgr->begin(), end = mdcMgr->end(); it != end;
3788 if ( it->helix[2] ) rPt[index] = 1 / fabs( it->helix[2] );
3789 if (
debug_ == 4 ) cout <<
"rPt...[ " << index <<
" ]...." << rPt[index] << endl;
3790 if ( rPt[index] < 0 ) rPt[index] =
DBL_MAX;
3792 std::vector<MdcRec_wirhit*> pt = it->hitcol;
3793 if (
debug_ == 4 ) cout <<
"ppt size: " << pt.size() << endl;
3794 int outermost( -1 );
3795 for ( vector<MdcRec_wirhit*>::iterator ii = pt.end() - 1; ii != pt.begin() - 1; ii-- )
3797 int lyr( ( *ii )->geo->Lyr()->Id() );
3798 if ( outermost < lyr ) outermost = lyr;
3799 if (
debug_ == 4 ) cout <<
"outmost: " << outermost <<
" lyr: " << lyr << endl;
3801 rOM[index] = outermost;
3802 order[index] = index;
3807 for (
int j, k = ntrk - 1; k >= 0; k = j )
3810 for (
int i = 1; i <= k; i++ )
3811 if ( rPt[i - 1] < rPt[i] )
3814 std::swap( order[i], order[j] );
3815 std::swap( rPt[i], rPt[j] );
3816 std::swap( rOM[i], rOM[j] );
3817 std::swap( rCont[i], rCont[j] );
3818 std::swap( rGen[i], rGen[j] );
3825 DataObject* aReconEvent;
3826 eventSvc()->findObject(
"/Event/Recon", aReconEvent );
3831 StatusCode sc = eventSvc()->registerObject(
"/Event/Recon", recevt );
3832 if ( sc != StatusCode::SUCCESS )
3834 log << MSG::FATAL <<
"Could not register ReconEvent" << endmsg;
3842 log << MSG::INFO <<
"beginning to make RecMdcKalTrackCol" << endmsg;
3845 for (
int l = 0; l < ntrk; l++ )
3848 MdcRec_trk& TrasanTRK = *( mdcMgr->begin() + order[l] );
3849 MdcRec_trk_add& TrasanTRK_add = *( mdc_addMgr->begin() + order[l] );
3852 int trasqual = TrasanTRK_add.
quality;
3853 if (
debug_ == 4 ) cout <<
"kalman_fitting>trasqual... " << trasqual << endl;
3854 if ( trasqual )
continue;
3857 if (
debug_ == 4 ) cout <<
"******* KalFit NUMBER : " << newcount << std::endl;
3861 if ( ( TrasanTRK_add.
decision & 32 ) == 32 || ( TrasanTRK_add.
decision & 64 ) == 64 )
3868 for (
int i = 0; i < 5; i++ ) a[i] = TrasanTRK.
helix[i];
3870 HepSymMatrix ea( 5 );
3871 for (
int i = 0, k = 0; i < 5; i++ )
3873 for (
int j = 0; j <= i; j++ )
3876 ea[j][i] = ea[i][j];
3882 double fiTerm = TrasanTRK.
fiTerm;
3888 int inlyr( 999 ), outlyr( -1 );
3889 int* rStat =
new int[43];
3890 for (
int irStat = 0; irStat < 43; ++irStat ) rStat[irStat] = 0;
3891 std::vector<MdcRec_wirhit*> pt = TrasanTRK.
hitcol;
3893 if (
debug_ == 4 ) cout <<
"*********Pt size****" << pt.size() << endl;
3895 int Num[43] = { 0 };
3896 for ( vector<MdcRec_wirhit*>::iterator ii = pt.end() - 1; ii != pt.begin() - 1; ii-- )
3897 { Num[( *ii )->geo->Lyr()->Id()]++; }
3900 for ( vector<MdcRec_wirhit*>::iterator ii = pt.end() - 1; ii != pt.begin() - 1; ii-- )
3904 if ( Num[( *ii )->geo->Lyr()->Id()] > 3 )
3907 cout <<
"WARNING: I found " << Num[( *ii )->geo->Lyr()->Id()]
3908 <<
" hits in the layer " << ( *ii )->geo->Lyr()->Id() << std::endl;
3935 double dist[2] = { rechit.
ddl, rechit.
ddr };
3936 double erdist[2] = { rechit.
erddl, rechit.
erddr };
3939 int lr_decision( 0 );
3942 if ( rechit.
lr == 2 || rechit.
lr == 0 ) lr_decision = -1;
3944 else if ( rechit.
lr == 1 ) lr_decision = 1;
3947 int ind( geo->
Lyr()->
Id() );
3952 if ( inlyr > ind ) inlyr = ind;
3953 if ( outlyr < ind ) outlyr = ind;
3955 if (
debug_ == 4 ) cout <<
"**** NUMBER OF Mdc HITS (TRASAN) = " << hit_asso << std::endl;
3958 int empty_between( 0 ), empty( 0 );
3959 for (
int i = inlyr; i <= outlyr; i++ )
3960 if ( !rStat[i] ) empty_between++;
3961 empty = empty_between + inlyr + ( 42 - outlyr );
3967 unsigned int nhit = track_lead.
HitsMdc().size();
3968 if ( !nhit &&
debug_ == 4 )
3970 cout <<
" ATTENTION TRACK WITH ONLY HITS " << nhit << std::endl;
3975 double KalFitst( 0 ), KalFitax( 0 ), KalFitschi2( 0 );
3977 Hep3Vector outer_pivot( track_lead.
x( fiTerm ) );
3981 std::cout <<
"before track_lead.pivot(outer_pivot) ,the error matrix of track_lead is .."
3982 << track_lead.
Ea() << std::endl;
3984 track_lead.
pivot( outer_pivot );
3988 HepSymMatrix Eakal( 5, 0 );
3992 double costheta = track_lead.
a()[4] / sqrt( 1.0 + track_lead.
a()[4] * track_lead.
a()[4] );
3993 if ( ( 1.0 / fabs( track_lead.
a()[2] ) <
pt_cut_ ) && ( fabs( costheta ) >
theta_cut_ ) )
4002 std::cout <<
"from Mdc Pattern Recognition: " << std::endl;
4004 Helix work( track_lead.
pivot(), track_lead.
a(), track_lead.
Ea() );
4006 std::cout <<
" dr = " << work.
a()[0] <<
", Er_dr = " << sqrt( work.
Ea()[0][0] )
4008 std::cout <<
" phi0 = " << work.
a()[1] <<
", Er_phi0 = " << sqrt( work.
Ea()[1][1] )
4010 std::cout <<
" PT = " << 1 / work.
a()[2] <<
", Er_kappa = " << sqrt( work.
Ea()[2][2] )
4012 std::cout <<
" dz = " << work.
a()[3] <<
", Er_dz = " << sqrt( work.
Ea()[3][3] )
4014 std::cout <<
" tanl = " << work.
a()[4] <<
", Er_tanl = " << sqrt( work.
Ea()[4][4] )
4024 cout <<
" Mdc FIRST KALMAN FIT " << std::endl;
4025 Helix work( track_lead.
pivot(), track_lead.
a(), track_lead.
Ea() );
4027 cout <<
" dr = " << work.
a()[0] <<
", Er_dr = " << sqrt( work.
Ea()[0][0] ) << std::endl;
4028 cout <<
" phi0 = " << work.
a()[1] <<
", Er_phi0 = " << sqrt( work.
Ea()[1][1] )
4030 cout <<
" PT = " << 1 / work.
a()[2] <<
", Er_kappa = " << sqrt( work.
Ea()[2][2] )
4032 cout <<
" dz = " << work.
a()[3] <<
", Er_dz = " << sqrt( work.
Ea()[3][3] ) << std::endl;
4033 cout <<
" tanl = " << work.
a()[4] <<
", Er_tanl = " << sqrt( work.
Ea()[4][4] )
4041 complete_track( TrasanTRK, TrasanTRK_add, track_lead, kaltrk, kalcol, segcol );
4046 DataObject* aRecKalEvent;
4047 eventSvc()->findObject(
"/Event/Recon/RecMdcKalTrackCol", aRecKalEvent );
4051 kalsc = eventSvc()->unregisterObject(
"/Event/Recon/RecMdcKalTrackCol" );
4052 if ( kalsc != StatusCode::SUCCESS )
4054 log << MSG::FATAL <<
"Could not unregister RecMdcKalTrack collection" << endmsg;
4059 kalsc = eventSvc()->registerObject(
"/Event/Recon/RecMdcKalTrackCol", kalcol );
4060 if ( kalsc.isFailure() )
4062 log << MSG::FATAL <<
"Could not register RecMdcKalTrack" << endmsg;
4065 log << MSG::INFO <<
"RecMdcKalTrackCol registered successfully!" << endmsg;
4069 DataObject* aRecKalSegEvent;
4070 eventSvc()->findObject(
"/Event/Recon/RecMdcKalHelixSegCol", aRecKalSegEvent );
4071 if ( aRecKalSegEvent )
4074 segsc = eventSvc()->unregisterObject(
"/Event/Recon/RecMdcKalHelixSegCol" );
4075 if ( segsc != StatusCode::SUCCESS )
4077 log << MSG::FATAL <<
"Could not unregister RecMdcKalHelixSegCol collection" << endmsg;
4082 segsc = eventSvc()->registerObject(
"/Event/Recon/RecMdcKalHelixSegCol", segcol );
4083 if ( segsc.isFailure() )
4085 log << MSG::FATAL <<
"Could not register RecMdcKalHelixSeg" << endmsg;
4088 log << MSG::INFO <<
"RecMdcKalHelixSegCol registered successfully!" << endmsg;
4090 double x1( 0. ), x2( 0. ), y1( 0. ), y2( 0. ), z1( 0. ), z2( 0. ), the1( 0. ), the2( 0. ),
4092 double r1( 0. ), r2( 0. ), kap1( 999. ), kap2( 999. ), tanl1( 0. ), tanl2( 0. );
4095 SmartDataPtr<RecMdcKalTrackCol> kaltrkCol( eventSvc(),
"/Event/Recon/RecMdcKalTrackCol" );
4098 log << MSG::FATAL <<
"Could not find RecMdcKalTrackCol" << endmsg;
4101 log << MSG::INFO <<
"Begin to check RecMdcKalTrackCol" << endmsg;
4102 RecMdcKalTrackCol::iterator iter_trk = kaltrkCol->begin();
4103 for (
int jj = 1; iter_trk != kaltrkCol->end(); iter_trk++, jj++ )
4105 log << MSG::DEBUG <<
"retrieved MDC Kalmantrack:"
4106 <<
"Track Id: " << ( *iter_trk )->getTrackId()
4107 <<
" Mass of the fit: " << ( *iter_trk )->getMass( 2 ) << endmsg
4108 <<
" Length of the track: " << ( *iter_trk )->getLength( 2 )
4109 <<
" Tof of the track: " << ( *iter_trk )->getTof( 2 ) << endmsg
4110 <<
" Chisq of the fit: " << ( *iter_trk )->getChisq( 0, 2 ) <<
" "
4111 << ( *iter_trk )->getChisq( 1, 2 ) << endmsg
4112 <<
"Ndf of the fit: " << ( *iter_trk )->getNdf( 0, 1 ) <<
" "
4113 << ( *iter_trk )->getNdf( 1, 2 ) << endmsg <<
"Kappa " << ( *iter_trk )->getZHelix()[2]
4118 { std::cout <<
"the size of gothelixsegs ..." << gothelixsegs.size() << std::endl; }
4120 HelixSegRefVec::iterator it_gothelixseg = gothelixsegs.begin();
4121 for ( ; it_gothelixseg != gothelixsegs.end(); it_gothelixseg++ )
4125 std::cout <<
"the layerId of this helixseg is ..." << ( *it_gothelixseg )->getLayerId()
4127 std::cout <<
"the residual of this helixseg exclude the meas hit"
4128 << ( *it_gothelixseg )->getResExcl() << std::endl;
4129 std::cout <<
"the residual of this helixseg include the meas hit"
4130 << ( *it_gothelixseg )->getResIncl() << std::endl;
4131 std::cout <<
"the track id of the helixseg is ..." << ( *it_gothelixseg )->getTrackId()
4133 std::cout <<
"the tof of the helixseg is ..." << ( *it_gothelixseg )->getTof()
4135 std::cout <<
"the Zhit of the helixseg is ..." << ( *it_gothelixseg )->getZhit()
4139 for (
int i = 0; i < 43; i++ )
4141 log << MSG::DEBUG <<
"retrieved pathl[" << i <<
"]= " << ( *iter_trk )->getPathl( i )
4147 m_trackid = ( *iter_trk )->getTrackId();
4148 for (
int jj = 0, iii = 0; jj < 5; jj++ )
4150 m_length[jj] = ( *iter_trk )->getLength( jj );
4151 m_tof[jj] = ( *iter_trk )->getTof( jj );
4152 m_nhits[jj] = ( *iter_trk )->getNhits( jj );
4153 m_zhelix[jj] = ( *iter_trk )->getZHelix()[jj];
4154 m_zhelixe[jj] = ( *iter_trk )->getZHelixE()[jj];
4155 m_zhelixmu[jj] = ( *iter_trk )->getZHelixMu()[jj];
4156 m_zhelixk[jj] = ( *iter_trk )->getZHelixK()[jj];
4157 m_zhelixp[jj] = ( *iter_trk )->getZHelixP()[jj];
4158 m_fhelix[jj] = ( *iter_trk )->getFHelix()[jj];
4159 m_fhelixe[jj] = ( *iter_trk )->getFHelixE()[jj];
4160 m_fhelixmu[jj] = ( *iter_trk )->getFHelixMu()[jj];
4161 m_fhelixk[jj] = ( *iter_trk )->getFHelixK()[jj];
4162 m_fhelixp[jj] = ( *iter_trk )->getFHelixP()[jj];
4163 m_lhelix[jj] = ( *iter_trk )->getLHelix()[jj];
4164 m_lhelixe[jj] = ( *iter_trk )->getLHelixE()[jj];
4165 m_lhelixmu[jj] = ( *iter_trk )->getLHelixMu()[jj];
4166 m_lhelixk[jj] = ( *iter_trk )->getLHelixK()[jj];
4167 m_lhelixp[jj] = ( *iter_trk )->getLHelixP()[jj];
4170 for (
int kk = 0; kk <= jj; kk++, iii++ )
4172 m_zerror[iii] = ( *iter_trk )->getZError()[jj][kk];
4173 m_zerrore[iii] = ( *iter_trk )->getZErrorE()[jj][kk];
4174 m_zerrormu[iii] = ( *iter_trk )->getZErrorMu()[jj][kk];
4175 m_zerrork[iii] = ( *iter_trk )->getZErrorK()[jj][kk];
4176 m_zerrorp[iii] = ( *iter_trk )->getZErrorP()[jj][kk];
4177 m_ferror[iii] = ( *iter_trk )->getFError()[jj][kk];
4178 m_ferrore[iii] = ( *iter_trk )->getFErrorE()[jj][kk];
4179 m_ferrormu[iii] = ( *iter_trk )->getFErrorMu()[jj][kk];
4180 m_ferrork[iii] = ( *iter_trk )->getFErrorK()[jj][kk];
4181 m_ferrorp[iii] = ( *iter_trk )->getFErrorP()[jj][kk];
4182 m_lerror[iii] = ( *iter_trk )->getLError()[jj][kk];
4183 m_lerrore[iii] = ( *iter_trk )->getLErrorE()[jj][kk];
4184 m_lerrormu[iii] = ( *iter_trk )->getLErrorMu()[jj][kk];
4185 m_lerrork[iii] = ( *iter_trk )->getLErrorK()[jj][kk];
4186 m_lerrorp[iii] = ( *iter_trk )->getLErrorP()[jj][kk];
4225 m_chisq[0][0] = ( *iter_trk )->getChisq( 0, 0 );
4226 m_chisq[1][0] = ( *iter_trk )->getChisq( 0, 1 );
4227 m_chisq[2][0] = ( *iter_trk )->getChisq( 0, 2 );
4228 m_chisq[3][0] = ( *iter_trk )->getChisq( 0, 3 );
4229 m_chisq[4][0] = ( *iter_trk )->getChisq( 0, 4 );
4230 m_chisq[0][1] = ( *iter_trk )->getChisq( 1, 0 );
4231 m_chisq[1][1] = ( *iter_trk )->getChisq( 1, 1 );
4232 m_chisq[2][1] = ( *iter_trk )->getChisq( 1, 2 );
4233 m_chisq[3][1] = ( *iter_trk )->getChisq( 1, 3 );
4234 m_chisq[4][1] = ( *iter_trk )->getChisq( 1, 4 );
4236 m_ndf[0][0] = ( *iter_trk )->getNdf( 0, 0 );
4237 m_ndf[1][0] = ( *iter_trk )->getNdf( 0, 1 );
4238 m_ndf[2][0] = ( *iter_trk )->getNdf( 0, 2 );
4239 m_ndf[3][0] = ( *iter_trk )->getNdf( 0, 3 );
4240 m_ndf[4][0] = ( *iter_trk )->getNdf( 0, 4 );
4241 m_ndf[0][1] = ( *iter_trk )->getNdf( 1, 0 );
4242 m_ndf[1][1] = ( *iter_trk )->getNdf( 1, 1 );
4243 m_ndf[2][1] = ( *iter_trk )->getNdf( 1, 2 );
4244 m_ndf[3][1] = ( *iter_trk )->getNdf( 1, 3 );
4245 m_ndf[4][1] = ( *iter_trk )->getNdf( 1, 4 );
4247 m_stat[0][0] = ( *iter_trk )->getStat( 0, 0 );
4248 m_stat[1][0] = ( *iter_trk )->getStat( 0, 1 );
4249 m_stat[2][0] = ( *iter_trk )->getStat( 0, 2 );
4250 m_stat[3][0] = ( *iter_trk )->getStat( 0, 3 );
4251 m_stat[4][0] = ( *iter_trk )->getStat( 0, 4 );
4252 m_stat[0][1] = ( *iter_trk )->getStat( 1, 0 );
4253 m_stat[1][1] = ( *iter_trk )->getStat( 1, 1 );
4254 m_stat[2][1] = ( *iter_trk )->getStat( 1, 2 );
4255 m_stat[3][1] = ( *iter_trk )->getStat( 1, 3 );
4256 m_stat[4][1] = ( *iter_trk )->getStat( 1, 4 );
4258 m_fptot = sqrt( 1 + pow( m_fhelix[4], 2 ) ) / m_fhelix[2];
4259 m_fptote = sqrt( 1 + pow( m_fhelixe[4], 2 ) ) / m_fhelixe[2];
4260 m_fptotmu = sqrt( 1 + pow( m_fhelixmu[4], 2 ) ) / m_fhelixmu[2];
4261 m_fptotk = sqrt( 1 + pow( m_fhelixk[4], 2 ) ) / m_fhelixk[2];
4262 m_fptotp = sqrt( 1 + pow( m_fhelixp[4], 2 ) ) / m_fhelixp[2];
4264 m_zpt = 1 / m_zhelix[2];
4265 m_zpte = 1 / m_zhelixe[2];
4266 m_zptmu = 1 / m_zhelixmu[2];
4267 m_zptk = 1 / m_zhelixk[2];
4268 m_zptp = 1 / m_zhelixp[2];
4270 m_fpt = 1 / m_fhelix[2];
4271 m_fpte = 1 / m_fhelixe[2];
4272 m_fptmu = 1 / m_fhelixmu[2];
4273 m_fptk = 1 / m_fhelixk[2];
4274 m_fptp = 1 / m_fhelixp[2];
4276 m_lpt = 1 / m_lhelix[2];
4277 m_lpte = 1 / m_lhelixe[2];
4278 m_lptmu = 1 / m_lhelixmu[2];
4279 m_lptk = 1 / m_lhelixk[2];
4280 m_lptp = 1 / m_lhelixp[2];
4282 m_lptot = sqrt( 1 + pow( m_lhelix[4], 2 ) ) / m_lhelix[2];
4283 m_lptote = sqrt( 1 + pow( m_lhelixe[4], 2 ) ) / m_lhelixe[2];
4284 m_lptotmu = sqrt( 1 + pow( m_lhelixmu[4], 2 ) ) / m_lhelixmu[2];
4285 m_lptotk = sqrt( 1 + pow( m_lhelixk[4], 2 ) ) / m_lhelixk[2];
4286 m_lptotp = sqrt( 1 + pow( m_lhelixp[4], 2 ) ) / m_lhelixp[2];
4288 m_zptot = sqrt( 1 + pow( m_zhelix[4], 2 ) ) / m_zhelix[2];
4289 m_zptote = sqrt( 1 + pow( m_zhelixe[4], 2 ) ) / m_zhelixe[2];
4290 m_zptotmu = sqrt( 1 + pow( m_zhelixmu[4], 2 ) ) / m_zhelixmu[2];
4291 m_zptotk = sqrt( 1 + pow( m_zhelixk[4], 2 ) ) / m_zhelixk[2];
4292 m_zptotp = sqrt( 1 + pow( m_zhelixp[4], 2 ) ) / m_zhelixp[2];
4295 m_zsigp = sqrt( pow( ( m_zptot / m_zhelix[2] ), 2 ) * m_zerror[5] +
4296 pow( ( m_zhelix[4] / m_zptot ), 2 ) * pow( ( 1 / m_zhelix[2] ), 4 ) *
4298 2 * m_zhelix[4] * m_zerror[12] * pow( ( 1 / m_zhelix[2] ), 3 ) );
4299 m_zsigpe = sqrt( pow( ( m_zptote / m_zhelixe[2] ), 2 ) * m_zerrore[5] +
4300 pow( ( m_zhelixe[4] / m_zptote ), 2 ) *
4301 pow( ( 1 / m_zhelixe[2] ), 4 ) * m_zerrore[14] -
4302 2 * m_zhelixe[4] * m_zerrore[12] * pow( ( 1 / m_zhelixe[2] ), 3 ) );
4304 sqrt( pow( ( m_zptotmu / m_zhelixmu[2] ), 2 ) * m_zerrormu[5] +
4305 pow( ( m_zhelixmu[4] / m_zptotmu ), 2 ) * pow( ( 1 / m_zhelixmu[2] ), 4 ) *
4307 2 * m_zhelixmu[4] * m_zerrormu[12] * pow( ( 1 / m_zhelixmu[2] ), 3 ) );
4308 m_zsigpk = sqrt( pow( ( m_zptotk / m_zhelixk[2] ), 2 ) * m_zerrork[5] +
4309 pow( ( m_zhelixk[4] / m_zptotk ), 2 ) *
4310 pow( ( 1 / m_zhelixk[2] ), 4 ) * m_zerrork[14] -
4311 2 * m_zhelixk[4] * m_zerrork[12] * pow( ( 1 / m_zhelixk[2] ), 3 ) );
4312 m_zsigpp = sqrt( pow( ( m_zptotp / m_zhelixp[2] ), 2 ) * m_zerrorp[5] +
4313 pow( ( m_zhelixp[4] / m_zptotp ), 2 ) *
4314 pow( ( 1 / m_zhelixp[2] ), 4 ) * m_zerrorp[14] -
4315 2 * m_zhelixp[4] * m_zerrorp[12] * pow( ( 1 / m_zhelixp[2] ), 3 ) );
4318 StatusCode sc1 = m_nt1->write();
4319 if ( sc1.isFailure() ) cout <<
"Ntuple1 filling failed!" << endl;
4326 phi1 = ( *iter_trk )->getFFi0();
4327 r1 = ( *iter_trk )->getFDr();
4328 z1 = ( *iter_trk )->getFDz();
4329 kap1 = ( *iter_trk )->getFCpa();
4330 tanl1 = ( *iter_trk )->getFTanl();
4333 p1 = sqrt( 1 + tanl1 * tanl1 ) / kap1;
4334 the1 =
M_PI / 2 - atan( tanl1 );
4338 phi2 = ( *iter_trk )->getFFi0();
4339 r2 = ( *iter_trk )->getFDr();
4340 z2 = ( *iter_trk )->getFDz();
4341 kap2 = ( *iter_trk )->getFCpa();
4342 tanl2 = ( *iter_trk )->getFTanl();
4345 p2 = sqrt( 1 + tanl2 * tanl2 ) / kap1;
4346 the2 =
M_PI / 2 - atan( tanl2 );
4355 m_delthe = the1 + the2;
4358 StatusCode sc2 = m_nt2->write();
4359 if ( sc2.isFailure() ) cout <<
"Ntuple2 filling failed!" << endl;
4366 if (
debug_ == 4 ) cout <<
"Kalfitting finished " << std::endl;
4372 MsgStream log(
msgSvc(), name() );
4373 double Pt_threshold( 0.3 );
4374 Hep3Vector IP( 0, 0, 0 );
4381 if ( !&whMgr )
return;
4384 int ntrk = mdcMgr->size();
4387 int nhits = whMgr->size();
4391 DataObject* aReconEvent;
4392 eventSvc()->findObject(
"/Event/Recon", aReconEvent );
4397 StatusCode sc = eventSvc()->registerObject(
"/Event/Recon", recevt );
4398 if ( sc != StatusCode::SUCCESS )
4400 log << MSG::FATAL <<
"Could not register ReconEvent" << endmsg;
4408 log << MSG::INFO <<
"beginning to make RecMdcKalTrackCol" << endmsg;
4410 MdcRec_trk& TrasanTRK = *( mdcMgr->begin() );
4419 if ( ( TrasanTRK_add.
decision & 32 ) == 32 || ( TrasanTRK_add.
decision & 64 ) == 64 )
4426 for (
int i = 0; i < 5; i++ ) a[i] = TrasanTRK.
helix[i];
4428 HepSymMatrix ea( 5 );
4429 for (
int i = 0, k = 0; i < 5; i++ )
4431 for (
int j = 0; j <= i; j++ )
4434 ea[j][i] = ea[i][j];
4440 double fiTerm = TrasanTRK.
fiTerm;
4448 int trasqual = TrasanTRK_add.
quality;
4449 if (
debug_ == 4 ) cout <<
"kalman_fitting>trasqual... " << trasqual << endl;
4450 if ( trasqual )
return;
4452 int inlyr( 999 ), outlyr( -1 );
4453 int* rStat =
new int[43];
4454 for (
int irStat = 0; irStat < 43; ++irStat ) rStat[irStat] = 0;
4455 std::vector<MdcRec_wirhit*> pt = TrasanTRK.
hitcol;
4457 if (
debug_ == 4 ) cout <<
"*********Pt size****" << pt.size() << endl;
4459 int Num[43] = { 0 };
4460 for ( vector<MdcRec_wirhit*>::iterator ii = pt.end() - 1; ii != pt.begin() - 1; ii-- )
4461 { Num[( *ii )->geo->Lyr()->Id()]++; }
4463 for ( vector<MdcRec_wirhit*>::iterator ii = pt.end() - 1; ii != pt.begin() - 1; ii-- )
4467 if ( Num[( *ii )->geo->Lyr()->Id()] > 3 )
4470 cout <<
"WARNING: I found " << Num[( *ii )->geo->Lyr()->Id()] <<
" hits in the layer "
4471 << ( *ii )->geo->Lyr()->Id() << std::endl;
4477 double dist[2] = { rechit.
ddl, rechit.
ddr };
4478 double erdist[2] = { rechit.
erddl, rechit.
erddr };
4481 int lr_decision( 0 );
4484 if ( rechit.
lr == 2 || rechit.
lr == 0 ) lr_decision = -1;
4486 else if ( rechit.
lr == 1 ) lr_decision = 1;
4489 int ind( geo->
Lyr()->
Id() );
4494 if ( inlyr > ind ) inlyr = ind;
4495 if ( outlyr < ind ) outlyr = ind;
4498 int empty_between( 0 ), empty( 0 );
4499 for (
int i = inlyr; i <= outlyr; i++ )
4500 if ( !rStat[i] ) empty_between++;
4501 empty = empty_between + inlyr + ( 42 - outlyr );
4504 if (
debug_ == 4 ) cout <<
"**** NUMBER OF Mdc HITS (TRASAN) = " << hit_asso << std::endl;
4509 unsigned int nhit = track_lead.
HitsMdc().size();
4512 cout <<
" ATTENTION TRACK WITH ONLY HITS " << nhit << std::endl;
4517 double KalFitst( 0 ), KalFitax( 0 ), KalFitschi2( 0 );
4519 Hep3Vector outer_pivot( track_lead.
x( fiTerm ) );
4523 std::cout <<
"before track_lead.pivot(outer_pivot) ,the error matrix of track_lead is .."
4524 << track_lead.
Ea() << std::endl;
4526 track_lead.
pivot( outer_pivot );
4530 HepSymMatrix Eakal( 5, 0 );
4534 double costheta = track_lead.
a()[4] / sqrt( 1.0 + track_lead.
a()[4] * track_lead.
a()[4] );
4535 if ( ( 1.0 / fabs( track_lead.
a()[2] ) <
pt_cut_ ) && ( fabs( costheta ) >
theta_cut_ ) )
4544 std::cout <<
"from Mdc Pattern Recognition: " << std::endl;
4546 Helix work( track_lead.
pivot(), track_lead.
a(), track_lead.
Ea() );
4548 std::cout <<
" dr = " << work.
a()[0] <<
", Er_dr = " << sqrt( work.
Ea()[0][0] )
4550 std::cout <<
" phi0 = " << work.
a()[1] <<
", Er_phi0 = " << sqrt( work.
Ea()[1][1] )
4552 std::cout <<
" PT = " << 1 / work.
a()[2] <<
", Er_kappa = " << sqrt( work.
Ea()[2][2] )
4554 std::cout <<
" dz = " << work.
a()[3] <<
", Er_dz = " << sqrt( work.
Ea()[3][3] )
4556 std::cout <<
" tanl = " << work.
a()[4] <<
", Er_tanl = " << sqrt( work.
Ea()[4][4] )
4565 cout <<
" Mdc FIRST KALMAN FIT " << std::endl;
4566 Helix work1( track_lead.
pivot(), track_lead.
a(), track_lead.
Ea() );
4568 cout <<
" dr = " << work1.
a()[0] <<
", Er_dr = " << sqrt( work1.
Ea()[0][0] ) << std::endl;
4569 cout <<
" phi0 = " << work1.
a()[1] <<
", Er_phi0 = " << sqrt( work1.
Ea()[1][1] )
4571 cout <<
" PT = " << 1 / work1.
a()[2] <<
", Er_kappa = " << sqrt( work1.
Ea()[2][2] )
4573 cout <<
" dz = " << work1.
a()[3] <<
", Er_dz = " << sqrt( work1.
Ea()[3][3] ) << std::endl;
4574 cout <<
" tanl = " << work1.
a()[4] <<
", Er_tanl = " << sqrt( work1.
Ea()[4][4] )
4582 complete_track( TrasanTRK, TrasanTRK_add, track_lead, kaltrk, kalcol, segcol );
4586 DataObject* aRecKalEvent;
4587 eventSvc()->findObject(
"/Event/Recon/RecMdcKalTrackCol", aRecKalEvent );
4591 kalsc = eventSvc()->unregisterObject(
"/Event/Recon/RecMdcKalTrackCol" );
4592 if ( kalsc != StatusCode::SUCCESS )
4594 log << MSG::FATAL <<
"Could not unregister RecMdcKalTrack collection" << endmsg;
4599 kalsc = eventSvc()->registerObject(
"/Event/Recon/RecMdcKalTrackCol", kalcol );
4600 if ( kalsc.isFailure() )
4602 log << MSG::FATAL <<
"Could not register RecMdcKalTrack" << endmsg;
4605 log << MSG::INFO <<
"RecMdcKalTrackCol registered successfully!" << endmsg;
4609 DataObject* aRecKalSegEvent;
4610 eventSvc()->findObject(
"/Event/Recon/RecMdcKalHelixSegCol", aRecKalSegEvent );
4611 if ( aRecKalSegEvent )
4614 segsc = eventSvc()->unregisterObject(
"/Event/Recon/RecMdcKalHelixSegCol" );
4615 if ( segsc != StatusCode::SUCCESS )
4617 log << MSG::FATAL <<
"Could not unregister RecMdcKalHelixSegCol collection" << endmsg;
4622 segsc = eventSvc()->registerObject(
"/Event/Recon/RecMdcKalHelixSegCol", segcol );
4623 if ( segsc.isFailure() )
4625 log << MSG::FATAL <<
"Could not register RecMdcKalHelixSeg" << endmsg;
4628 log << MSG::INFO <<
"RecMdcKalHelixSegCol registered successfully!" << endmsg;
4630 double x1( 0. ), x2( 0. ), y1( 0. ), y2( 0. ), z1( 0. ), z2( 0. ), the1( 0. ), the2( 0. ),
4632 double r1( 0. ), r2( 0. ), kap1( 999. ), kap2( 999. ), tanl1( 0. ), tanl2( 0. );
4635 SmartDataPtr<RecMdcKalTrackCol> kaltrkCol( eventSvc(),
"/Event/Recon/RecMdcKalTrackCol" );
4638 log << MSG::FATAL <<
"Could not find RecMdcKalTrackCol" << endmsg;
4641 log << MSG::INFO <<
"Begin to check RecMdcKalTrackCol" << endmsg;
4642 RecMdcKalTrackCol::iterator iter_trk = kaltrkCol->begin();
4643 for (
int jj = 1; iter_trk != kaltrkCol->end(); iter_trk++, jj++ )
4645 log << MSG::DEBUG <<
"retrieved MDC Kalmantrack:"
4646 <<
"Track Id: " << ( *iter_trk )->getTrackId()
4647 <<
" Mass of the fit: " << ( *iter_trk )->getMass( 2 ) << endmsg
4648 <<
" Length of the track: " << ( *iter_trk )->getLength( 2 )
4649 <<
" Tof of the track: " << ( *iter_trk )->getTof( 2 ) << endmsg
4650 <<
" Chisq of the fit: " << ( *iter_trk )->getChisq( 0, 2 ) <<
" "
4651 << ( *iter_trk )->getChisq( 1, 2 ) << endmsg
4652 <<
"Ndf of the fit: " << ( *iter_trk )->getNdf( 0, 1 ) <<
" "
4653 << ( *iter_trk )->getNdf( 1, 2 ) << endmsg <<
"Kappa " << ( *iter_trk )->getZHelix()[2]
4654 <<
"zhelixmu " << ( *iter_trk )->getZHelixMu() << endmsg;
4658 { std::cout <<
"the size of gothelixsegs ..." << gothelixsegs.size() << std::endl; }
4660 HelixSegRefVec::iterator it_gothelixseg = gothelixsegs.begin();
4661 for ( ; it_gothelixseg != gothelixsegs.end(); it_gothelixseg++ )
4665 std::cout <<
"the layerId of this helixseg is ..." << ( *it_gothelixseg )->getLayerId()
4667 std::cout <<
"the residual of this helixseg exclude the meas hit"
4668 << ( *it_gothelixseg )->getResExcl() << std::endl;
4669 std::cout <<
"the residual of this helixseg include the meas hit"
4670 << ( *it_gothelixseg )->getResIncl() << std::endl;
4671 std::cout <<
"the track id of the helixseg is ..." << ( *it_gothelixseg )->getTrackId()
4673 std::cout <<
"the tof of the helixseg is ..." << ( *it_gothelixseg )->getTof()
4675 std::cout <<
"the Zhit of the helixseg is ..." << ( *it_gothelixseg )->getZhit()
4679 for (
int i = 0; i < 43; i++ )
4681 log << MSG::DEBUG <<
"retrieved pathl[" << i <<
"]= " << ( *iter_trk )->getPathl( i )
4687 m_trackid = ( *iter_trk )->getTrackId();
4688 for (
int jj = 0, iii = 0; jj < 5; jj++ )
4690 m_length[jj] = ( *iter_trk )->getLength( jj );
4691 m_tof[jj] = ( *iter_trk )->getTof( jj );
4692 m_nhits[jj] = ( *iter_trk )->getNhits( jj );
4693 m_zhelix[jj] = ( *iter_trk )->getZHelix()[jj];
4694 m_zhelixe[jj] = ( *iter_trk )->getZHelixE()[jj];
4695 m_zhelixmu[jj] = ( *iter_trk )->getZHelixMu()[jj];
4696 m_zhelixk[jj] = ( *iter_trk )->getZHelixK()[jj];
4697 m_zhelixp[jj] = ( *iter_trk )->getZHelixP()[jj];
4698 m_fhelix[jj] = ( *iter_trk )->getFHelix()[jj];
4699 m_fhelixe[jj] = ( *iter_trk )->getFHelixE()[jj];
4700 m_fhelixmu[jj] = ( *iter_trk )->getFHelixMu()[jj];
4701 m_fhelixk[jj] = ( *iter_trk )->getFHelixK()[jj];
4702 m_fhelixp[jj] = ( *iter_trk )->getFHelixP()[jj];
4703 m_lhelix[jj] = ( *iter_trk )->getLHelix()[jj];
4704 m_lhelixe[jj] = ( *iter_trk )->getLHelixE()[jj];
4705 m_lhelixmu[jj] = ( *iter_trk )->getLHelixMu()[jj];
4706 m_lhelixk[jj] = ( *iter_trk )->getLHelixK()[jj];
4707 m_lhelixp[jj] = ( *iter_trk )->getLHelixP()[jj];
4710 for (
int kk = 0; kk <= jj; kk++, iii++ )
4712 m_zerror[iii] = ( *iter_trk )->getZError()[jj][kk];
4713 m_zerrore[iii] = ( *iter_trk )->getZErrorE()[jj][kk];
4714 m_zerrormu[iii] = ( *iter_trk )->getZErrorMu()[jj][kk];
4715 m_zerrork[iii] = ( *iter_trk )->getZErrorK()[jj][kk];
4716 m_zerrorp[iii] = ( *iter_trk )->getZErrorP()[jj][kk];
4717 m_ferror[iii] = ( *iter_trk )->getFError()[jj][kk];
4718 m_ferrore[iii] = ( *iter_trk )->getFErrorE()[jj][kk];
4719 m_ferrormu[iii] = ( *iter_trk )->getFErrorMu()[jj][kk];
4720 m_ferrork[iii] = ( *iter_trk )->getFErrorK()[jj][kk];
4721 m_ferrorp[iii] = ( *iter_trk )->getFErrorP()[jj][kk];
4722 m_lerror[iii] = ( *iter_trk )->getLError()[jj][kk];
4723 m_lerrore[iii] = ( *iter_trk )->getLErrorE()[jj][kk];
4724 m_lerrormu[iii] = ( *iter_trk )->getLErrorMu()[jj][kk];
4725 m_lerrork[iii] = ( *iter_trk )->getLErrorK()[jj][kk];
4726 m_lerrorp[iii] = ( *iter_trk )->getLErrorP()[jj][kk];
4732 m_chisq[0][0] = ( *iter_trk )->getChisq( 0, 0 );
4733 m_chisq[1][0] = ( *iter_trk )->getChisq( 0, 1 );
4734 m_chisq[2][0] = ( *iter_trk )->getChisq( 0, 2 );
4735 m_chisq[3][0] = ( *iter_trk )->getChisq( 0, 3 );
4736 m_chisq[4][0] = ( *iter_trk )->getChisq( 0, 4 );
4737 m_chisq[0][1] = ( *iter_trk )->getChisq( 1, 0 );
4738 m_chisq[1][1] = ( *iter_trk )->getChisq( 1, 1 );
4739 m_chisq[2][1] = ( *iter_trk )->getChisq( 1, 2 );
4740 m_chisq[3][1] = ( *iter_trk )->getChisq( 1, 3 );
4741 m_chisq[4][1] = ( *iter_trk )->getChisq( 1, 4 );
4743 m_ndf[0][0] = ( *iter_trk )->getNdf( 0, 0 );
4744 m_ndf[1][0] = ( *iter_trk )->getNdf( 0, 1 );
4745 m_ndf[2][0] = ( *iter_trk )->getNdf( 0, 2 );
4746 m_ndf[3][0] = ( *iter_trk )->getNdf( 0, 3 );
4747 m_ndf[4][0] = ( *iter_trk )->getNdf( 0, 4 );
4748 m_ndf[0][1] = ( *iter_trk )->getNdf( 1, 0 );
4749 m_ndf[1][1] = ( *iter_trk )->getNdf( 1, 1 );
4750 m_ndf[2][1] = ( *iter_trk )->getNdf( 1, 2 );
4751 m_ndf[3][1] = ( *iter_trk )->getNdf( 1, 3 );
4752 m_ndf[4][1] = ( *iter_trk )->getNdf( 1, 4 );
4754 m_stat[0][0] = ( *iter_trk )->getStat( 0, 0 );
4755 m_stat[1][0] = ( *iter_trk )->getStat( 0, 1 );
4756 m_stat[2][0] = ( *iter_trk )->getStat( 0, 2 );
4757 m_stat[3][0] = ( *iter_trk )->getStat( 0, 3 );
4758 m_stat[4][0] = ( *iter_trk )->getStat( 0, 4 );
4759 m_stat[0][1] = ( *iter_trk )->getStat( 1, 0 );
4760 m_stat[1][1] = ( *iter_trk )->getStat( 1, 1 );
4761 m_stat[2][1] = ( *iter_trk )->getStat( 1, 2 );
4762 m_stat[3][1] = ( *iter_trk )->getStat( 1, 3 );
4763 m_stat[4][1] = ( *iter_trk )->getStat( 1, 4 );
4765 m_fptot = sqrt( 1 + pow( m_fhelix[4], 2 ) ) / m_fhelix[2];
4766 m_fptote = sqrt( 1 + pow( m_fhelixe[4], 2 ) ) / m_fhelixe[2];
4767 m_fptotmu = sqrt( 1 + pow( m_fhelixmu[4], 2 ) ) / m_fhelixmu[2];
4768 m_fptotk = sqrt( 1 + pow( m_fhelixk[4], 2 ) ) / m_fhelixk[2];
4769 m_fptotp = sqrt( 1 + pow( m_fhelixp[4], 2 ) ) / m_fhelixp[2];
4771 m_zpt = 1 / m_zhelix[2];
4772 m_zpte = 1 / m_zhelixe[2];
4773 m_zptmu = 1 / m_zhelixmu[2];
4774 m_zptk = 1 / m_zhelixk[2];
4775 m_zptp = 1 / m_zhelixp[2];
4777 m_fpt = 1 / m_fhelix[2];
4778 m_fpte = 1 / m_fhelixe[2];
4779 m_fptmu = 1 / m_fhelixmu[2];
4780 m_fptk = 1 / m_fhelixk[2];
4781 m_fptp = 1 / m_fhelixp[2];
4783 m_lpt = 1 / m_lhelix[2];
4784 m_lpte = 1 / m_lhelixe[2];
4785 m_lptmu = 1 / m_lhelixmu[2];
4786 m_lptk = 1 / m_lhelixk[2];
4787 m_lptp = 1 / m_lhelixp[2];
4789 m_lptot = sqrt( 1 + pow( m_lhelix[4], 2 ) ) / m_lhelix[2];
4790 m_lptote = sqrt( 1 + pow( m_lhelixe[4], 2 ) ) / m_lhelixe[2];
4791 m_lptotmu = sqrt( 1 + pow( m_lhelixmu[4], 2 ) ) / m_lhelixmu[2];
4792 m_lptotk = sqrt( 1 + pow( m_lhelixk[4], 2 ) ) / m_lhelixk[2];
4793 m_lptotp = sqrt( 1 + pow( m_lhelixp[4], 2 ) ) / m_lhelixp[2];
4795 m_zptot = sqrt( 1 + pow( m_zhelix[4], 2 ) ) / m_zhelix[2];
4796 m_zptote = sqrt( 1 + pow( m_zhelixe[4], 2 ) ) / m_zhelixe[2];
4797 m_zptotmu = sqrt( 1 + pow( m_zhelixmu[4], 2 ) ) / m_zhelixmu[2];
4798 m_zptotk = sqrt( 1 + pow( m_zhelixk[4], 2 ) ) / m_zhelixk[2];
4799 m_zptotp = sqrt( 1 + pow( m_zhelixp[4], 2 ) ) / m_zhelixp[2];
4802 m_zsigp = sqrt( pow( ( m_zptot / m_zhelix[2] ), 2 ) * m_zerror[5] +
4803 pow( ( m_zhelix[4] / m_zptot ), 2 ) * pow( ( 1 / m_zhelix[2] ), 4 ) *
4805 2 * m_zhelix[4] * m_zerror[12] * pow( ( 1 / m_zhelix[2] ), 3 ) );
4806 m_zsigpe = sqrt( pow( ( m_zptote / m_zhelixe[2] ), 2 ) * m_zerrore[5] +
4807 pow( ( m_zhelixe[4] / m_zptote ), 2 ) *
4808 pow( ( 1 / m_zhelixe[2] ), 4 ) * m_zerrore[14] -
4809 2 * m_zhelixe[4] * m_zerrore[12] * pow( ( 1 / m_zhelixe[2] ), 3 ) );
4811 sqrt( pow( ( m_zptotmu / m_zhelixmu[2] ), 2 ) * m_zerrormu[5] +
4812 pow( ( m_zhelixmu[4] / m_zptotmu ), 2 ) * pow( ( 1 / m_zhelixmu[2] ), 4 ) *
4814 2 * m_zhelixmu[4] * m_zerrormu[12] * pow( ( 1 / m_zhelixmu[2] ), 3 ) );
4815 m_zsigpk = sqrt( pow( ( m_zptotk / m_zhelixk[2] ), 2 ) * m_zerrork[5] +
4816 pow( ( m_zhelixk[4] / m_zptotk ), 2 ) *
4817 pow( ( 1 / m_zhelixk[2] ), 4 ) * m_zerrork[14] -
4818 2 * m_zhelixk[4] * m_zerrork[12] * pow( ( 1 / m_zhelixk[2] ), 3 ) );
4819 m_zsigpp = sqrt( pow( ( m_zptotp / m_zhelixp[2] ), 2 ) * m_zerrorp[5] +
4820 pow( ( m_zhelixp[4] / m_zptotp ), 2 ) *
4821 pow( ( 1 / m_zhelixp[2] ), 4 ) * m_zerrorp[14] -
4822 2 * m_zhelixp[4] * m_zerrorp[12] * pow( ( 1 / m_zhelixp[2] ), 3 ) );
4825 StatusCode sc1 = m_nt1->write();
4826 if ( sc1.isFailure() ) cout <<
"Ntuple1 filling failed!" << endl;
4833 phi1 = ( *iter_trk )->getFFi0();
4834 r1 = ( *iter_trk )->getFDr();
4835 z1 = ( *iter_trk )->getFDz();
4836 kap1 = ( *iter_trk )->getFCpa();
4837 tanl1 = ( *iter_trk )->getFTanl();
4840 p1 = sqrt( 1 + tanl1 * tanl1 ) / kap1;
4841 the1 =
M_PI / 2 - atan( tanl1 );
4845 phi2 = ( *iter_trk )->getFFi0();
4846 r2 = ( *iter_trk )->getFDr();
4847 z2 = ( *iter_trk )->getFDz();
4848 kap2 = ( *iter_trk )->getFCpa();
4849 tanl2 = ( *iter_trk )->getFTanl();
4852 p2 = sqrt( 1 + tanl2 * tanl2 ) / kap1;
4853 the2 =
M_PI / 2 - atan( tanl2 );
4862 m_delthe = the1 + the2;
4865 StatusCode sc2 = m_nt2->write();
4866 if ( sc2.isFailure() ) cout <<
"Ntuple2 filling failed!" << endl;
4869 if (
debug_ == 4 ) cout <<
"Kalfitting finished " << std::endl;
4875 MsgStream log(
msgSvc(), name() );
4876 double Pt_threshold( 0.3 );
4877 Hep3Vector IP( 0, 0, 0 );
4884 if ( !&whMgr )
return;
4887 int ntrk = mdcMgr->size();
4890 int nhits = whMgr->size();
4893 double* rY =
new double[ntrk];
4894 double* rfiTerm =
new double[ntrk];
4895 double* rPt =
new double[ntrk];
4896 int* rOM =
new int[ntrk];
4897 unsigned int* order =
new unsigned int[ntrk];
4898 unsigned int* rCont =
new unsigned int[ntrk];
4899 unsigned int* rGen =
new unsigned int[ntrk];
4902 Hep3Vector csmp3[2];
4903 for ( vector<MdcRec_trk>::iterator it = mdcMgr->begin(), end = mdcMgr->end(); it != end;
4907 rfiTerm[index] = it->fiTerm;
4911 if ( it->helix[2] ) rPt[index] = 1 / fabs( it->helix[2] );
4912 if (
debug_ == 4 ) cout <<
"rPt...[ " << index <<
" ]...." << rPt[index] << endl;
4913 if ( rPt[index] < 0 ) rPt[index] =
DBL_MAX;
4915 std::vector<MdcRec_wirhit*> pt = it->hitcol;
4916 if (
debug_ == 4 ) cout <<
"ppt size: " << pt.size() << endl;
4917 int outermost( -1 );
4918 for ( vector<MdcRec_wirhit*>::iterator ii = pt.end() - 1; ii != pt.begin() - 1; ii-- )
4920 int lyr( ( *ii )->geo->Lyr()->Id() );
4921 if ( outermost < lyr )
4924 rY[index] = ( *ii )->geo->Forward().y();
4926 if (
debug_ == 4 ) cout <<
"outmost: " << outermost <<
" lyr: " << lyr << endl;
4928 rOM[index] = outermost;
4929 order[index] = index;
4934 for (
int j, k = ntrk - 1; k >= 0; k = j )
4937 for (
int i = 1; i <= k; i++ )
4938 if ( rY[i - 1] < rY[i] )
4941 std::swap( order[i], order[j] );
4942 std::swap( rY[i], rY[j] );
4943 std::swap( rOM[i], rOM[j] );
4944 std::swap( rCont[i], rCont[j] );
4945 std::swap( rGen[i], rGen[j] );
4954 DataObject* aReconEvent;
4955 eventSvc()->findObject(
"/Event/Recon", aReconEvent );
4960 StatusCode sc = eventSvc()->registerObject(
"/Event/Recon", recevt );
4961 if ( sc != StatusCode::SUCCESS )
4963 log << MSG::FATAL <<
"Could not register ReconEvent" << endmsg;
4971 log << MSG::INFO <<
"beginning to make RecMdcKalTrackCol" << endmsg;
4978 MdcRec_trk& TrasanTRK = *( mdcMgr->begin() + order[1] );
4979 MdcRec_trk_add& TrasanTRK_add = *( mdc_addMgr->begin() + order[1] );
4986 if (
debug_ == 4 ) cout <<
"******* KalFit NUMBER : " << newcount << std::endl;
4990 if ( ( TrasanTRK_add.
decision & 32 ) == 32 || ( TrasanTRK_add.
decision & 64 ) == 64 )
4997 for (
int i = 0; i < 5; i++ ) a[i] = TrasanTRK.
helix[i];
4999 HepSymMatrix ea( 5 );
5000 for (
int i = 0, k = 0; i < 5; i++ )
5002 for (
int j = 0; j <= i; j++ )
5005 ea[j][i] = ea[i][j];
5011 double fiTerm = TrasanTRK.
fiTerm;
5018 for (
int l = 0; l < ntrk; l++ )
5020 MdcRec_trk& TrasanTRK1 = *( mdcMgr->begin() + order[l] );
5021 MdcRec_trk_add& TrasanTRK_add1 = *( mdc_addMgr->begin() + order[l] );
5023 int trasqual = TrasanTRK_add1.
quality;
5024 if (
debug_ == 4 ) cout <<
"kalman_fitting>trasqual... " << trasqual << endl;
5025 if ( trasqual )
continue;
5027 int inlyr( 999 ), outlyr( -1 );
5028 int* rStat =
new int[43];
5029 for (
int irStat = 0; irStat < 43; ++irStat ) rStat[irStat] = 0;
5030 std::vector<MdcRec_wirhit*> pt = TrasanTRK1.
hitcol;
5032 if (
debug_ == 4 ) cout <<
"*********Pt size****" << pt.size() << endl;
5034 int Num[43] = { 0 };
5035 for ( vector<MdcRec_wirhit*>::iterator ii = pt.end() - 1; ii != pt.begin() - 1; ii-- )
5036 { Num[( *ii )->geo->Lyr()->Id()]++; }
5038 for ( vector<MdcRec_wirhit*>::iterator ii = pt.end() - 1; ii != pt.begin() - 1; ii-- )
5042 if ( Num[( *ii )->geo->Lyr()->Id()] > 3 )
5045 cout <<
"WARNING: I found " << Num[( *ii )->geo->Lyr()->Id()]
5046 <<
" hits in the layer " << ( *ii )->geo->Lyr()->Id() << std::endl;
5052 double dist[2] = { rechit.
ddl, rechit.
ddr };
5053 double erdist[2] = { rechit.
erddl, rechit.
erddr };
5056 int lr_decision( 0 );
5059 if ( rechit.
lr == 2 || rechit.
lr == 0 ) lr_decision = -1;
5061 else if ( rechit.
lr == 1 ) lr_decision = 1;
5064 int ind( geo->
Lyr()->
Id() );
5069 if ( inlyr > ind ) inlyr = ind;
5070 if ( outlyr < ind ) outlyr = ind;
5073 int empty_between( 0 ), empty( 0 );
5074 for (
int i = inlyr; i <= outlyr; i++ )
5075 if ( !rStat[i] ) empty_between++;
5076 empty = empty_between + inlyr + ( 42 - outlyr );
5079 if (
debug_ == 4 ) cout <<
"**** NUMBER OF Mdc HITS (TRASAN) = " << hit_asso << std::endl;
5084 unsigned int nhit = track_lead.
HitsMdc().size();
5087 cout <<
" ATTENTION TRACK WITH ONLY HITS " << nhit << std::endl;
5092 double KalFitst( 0 ), KalFitax( 0 ), KalFitschi2( 0 );
5094 Hep3Vector outer_pivot( track_lead.
x( fiTerm ) );
5098 std::cout <<
"before track_lead.pivot(outer_pivot) ,the error matrix of track_lead is .."
5099 << track_lead.
Ea() << std::endl;
5101 track_lead.
pivot( outer_pivot );
5105 HepSymMatrix Eakal( 5, 0 );
5109 double costheta = track_lead.
a()[4] / sqrt( 1.0 + track_lead.
a()[4] * track_lead.
a()[4] );
5110 if ( ( 1.0 / fabs( track_lead.
a()[2] ) <
pt_cut_ ) && ( fabs( costheta ) >
theta_cut_ ) )
5119 std::cout <<
"from Mdc Pattern Recognition: " << std::endl;
5121 Helix work( track_lead.
pivot(), track_lead.
a(), track_lead.
Ea() );
5123 std::cout <<
" dr = " << work.
a()[0] <<
", Er_dr = " << sqrt( work.
Ea()[0][0] )
5125 std::cout <<
" phi0 = " << work.
a()[1] <<
", Er_phi0 = " << sqrt( work.
Ea()[1][1] )
5127 std::cout <<
" PT = " << 1 / work.
a()[2] <<
", Er_kappa = " << sqrt( work.
Ea()[2][2] )
5129 std::cout <<
" dz = " << work.
a()[3] <<
", Er_dz = " << sqrt( work.
Ea()[3][3] )
5131 std::cout <<
" tanl = " << work.
a()[4] <<
", Er_tanl = " << sqrt( work.
Ea()[4][4] )
5140 cout <<
" Mdc FIRST KALMAN FIT " << std::endl;
5141 Helix work1( track_lead.
pivot(), track_lead.
a(), track_lead.
Ea() );
5143 cout <<
" dr = " << work1.
a()[0] <<
", Er_dr = " << sqrt( work1.
Ea()[0][0] ) << std::endl;
5144 cout <<
" phi0 = " << work1.
a()[1] <<
", Er_phi0 = " << sqrt( work1.
Ea()[1][1] )
5146 cout <<
" PT = " << 1 / work1.
a()[2] <<
", Er_kappa = " << sqrt( work1.
Ea()[2][2] )
5148 cout <<
" dz = " << work1.
a()[3] <<
", Er_dz = " << sqrt( work1.
Ea()[3][3] ) << std::endl;
5149 cout <<
" tanl = " << work1.
a()[4] <<
", Er_tanl = " << sqrt( work1.
Ea()[4][4] )
5157 complete_track( TrasanTRK, TrasanTRK_add, track_lead, kaltrk, kalcol, segcol );
5162 DataObject* aRecKalEvent;
5163 eventSvc()->findObject(
"/Event/Recon/RecMdcKalTrackCol", aRecKalEvent );
5167 kalsc = eventSvc()->unregisterObject(
"/Event/Recon/RecMdcKalTrackCol" );
5168 if ( kalsc != StatusCode::SUCCESS )
5170 log << MSG::FATAL <<
"Could not unregister RecMdcKalTrack collection" << endmsg;
5175 kalsc = eventSvc()->registerObject(
"/Event/Recon/RecMdcKalTrackCol", kalcol );
5176 if ( kalsc.isFailure() )
5178 log << MSG::FATAL <<
"Could not register RecMdcKalTrack" << endmsg;
5181 log << MSG::INFO <<
"RecMdcKalTrackCol registered successfully!" << endmsg;
5185 DataObject* aRecKalSegEvent;
5186 eventSvc()->findObject(
"/Event/Recon/RecMdcKalHelixSegCol", aRecKalSegEvent );
5187 if ( aRecKalSegEvent )
5190 segsc = eventSvc()->unregisterObject(
"/Event/Recon/RecMdcKalHelixSegCol" );
5191 if ( segsc != StatusCode::SUCCESS )
5193 log << MSG::FATAL <<
"Could not unregister RecMdcKalHelixSegCol collection" << endmsg;
5198 segsc = eventSvc()->registerObject(
"/Event/Recon/RecMdcKalHelixSegCol", segcol );
5199 if ( segsc.isFailure() )
5201 log << MSG::FATAL <<
"Could not register RecMdcKalHelixSeg" << endmsg;
5204 log << MSG::INFO <<
"RecMdcKalHelixSegCol registered successfully!" << endmsg;
5206 double x1( 0. ), x2( 0. ), y1( 0. ), y2( 0. ), z1( 0. ), z2( 0. ), the1( 0. ), the2( 0. ),
5208 double r1( 0. ), r2( 0. ), kap1( 999. ), kap2( 999. ), tanl1( 0. ), tanl2( 0. );
5211 SmartDataPtr<RecMdcKalTrackCol> kaltrkCol( eventSvc(),
"/Event/Recon/RecMdcKalTrackCol" );
5214 log << MSG::FATAL <<
"Could not find RecMdcKalTrackCol" << endmsg;
5217 log << MSG::INFO <<
"Begin to check RecMdcKalTrackCol" << endmsg;
5218 RecMdcKalTrackCol::iterator iter_trk = kaltrkCol->begin();
5219 for (
int jj = 1; iter_trk != kaltrkCol->end(); iter_trk++, jj++ )
5221 log << MSG::DEBUG <<
"retrieved MDC Kalmantrack:"
5222 <<
"Track Id: " << ( *iter_trk )->getTrackId()
5223 <<
" Mass of the fit: " << ( *iter_trk )->getMass( 2 ) << endmsg
5224 <<
" Length of the track: " << ( *iter_trk )->getLength( 2 )
5225 <<
" Tof of the track: " << ( *iter_trk )->getTof( 2 ) << endmsg
5226 <<
" Chisq of the fit: " << ( *iter_trk )->getChisq( 0, 2 ) <<
" "
5227 << ( *iter_trk )->getChisq( 1, 2 ) << endmsg
5228 <<
"Ndf of the fit: " << ( *iter_trk )->getNdf( 0, 1 ) <<
" "
5229 << ( *iter_trk )->getNdf( 1, 2 ) << endmsg <<
"Kappa " << ( *iter_trk )->getZHelix()[2]
5230 <<
"zhelixmu " << ( *iter_trk )->getZHelixMu() << endmsg;
5234 { std::cout <<
"the size of gothelixsegs ..." << gothelixsegs.size() << std::endl; }
5236 HelixSegRefVec::iterator it_gothelixseg = gothelixsegs.begin();
5237 for ( ; it_gothelixseg != gothelixsegs.end(); it_gothelixseg++ )
5241 std::cout <<
"the layerId of this helixseg is ..." << ( *it_gothelixseg )->getLayerId()
5243 std::cout <<
"the residual of this helixseg exclude the meas hit"
5244 << ( *it_gothelixseg )->getResExcl() << std::endl;
5245 std::cout <<
"the residual of this helixseg include the meas hit"
5246 << ( *it_gothelixseg )->getResIncl() << std::endl;
5247 std::cout <<
"the track id of the helixseg is ..." << ( *it_gothelixseg )->getTrackId()
5249 std::cout <<
"the tof of the helixseg is ..." << ( *it_gothelixseg )->getTof()
5251 std::cout <<
"the Zhit of the helixseg is ..." << ( *it_gothelixseg )->getZhit()
5255 for (
int i = 0; i < 43; i++ )
5257 log << MSG::DEBUG <<
"retrieved pathl[" << i <<
"]= " << ( *iter_trk )->getPathl( i )
5263 m_trackid = ( *iter_trk )->getTrackId();
5264 for (
int jj = 0, iii = 0; jj < 5; jj++ )
5266 m_length[jj] = ( *iter_trk )->getLength( jj );
5267 m_tof[jj] = ( *iter_trk )->getTof( jj );
5268 m_nhits[jj] = ( *iter_trk )->getNhits( jj );
5269 m_zhelix[jj] = ( *iter_trk )->getZHelix()[jj];
5270 m_zhelixe[jj] = ( *iter_trk )->getZHelixE()[jj];
5271 m_zhelixmu[jj] = ( *iter_trk )->getZHelixMu()[jj];
5272 m_zhelixk[jj] = ( *iter_trk )->getZHelixK()[jj];
5273 m_zhelixp[jj] = ( *iter_trk )->getZHelixP()[jj];
5274 m_fhelix[jj] = ( *iter_trk )->getFHelix()[jj];
5275 m_fhelixe[jj] = ( *iter_trk )->getFHelixE()[jj];
5276 m_fhelixmu[jj] = ( *iter_trk )->getFHelixMu()[jj];
5277 m_fhelixk[jj] = ( *iter_trk )->getFHelixK()[jj];
5278 m_fhelixp[jj] = ( *iter_trk )->getFHelixP()[jj];
5279 m_lhelix[jj] = ( *iter_trk )->getLHelix()[jj];
5280 m_lhelixe[jj] = ( *iter_trk )->getLHelixE()[jj];
5281 m_lhelixmu[jj] = ( *iter_trk )->getLHelixMu()[jj];
5282 m_lhelixk[jj] = ( *iter_trk )->getLHelixK()[jj];
5283 m_lhelixp[jj] = ( *iter_trk )->getLHelixP()[jj];
5286 for (
int kk = 0; kk <= jj; kk++, iii++ )
5288 m_zerror[iii] = ( *iter_trk )->getZError()[jj][kk];
5289 m_zerrore[iii] = ( *iter_trk )->getZErrorE()[jj][kk];
5290 m_zerrormu[iii] = ( *iter_trk )->getZErrorMu()[jj][kk];
5291 m_zerrork[iii] = ( *iter_trk )->getZErrorK()[jj][kk];
5292 m_zerrorp[iii] = ( *iter_trk )->getZErrorP()[jj][kk];
5293 m_ferror[iii] = ( *iter_trk )->getFError()[jj][kk];
5294 m_ferrore[iii] = ( *iter_trk )->getFErrorE()[jj][kk];
5295 m_ferrormu[iii] = ( *iter_trk )->getFErrorMu()[jj][kk];
5296 m_ferrork[iii] = ( *iter_trk )->getFErrorK()[jj][kk];
5297 m_ferrorp[iii] = ( *iter_trk )->getFErrorP()[jj][kk];
5298 m_lerror[iii] = ( *iter_trk )->getLError()[jj][kk];
5299 m_lerrore[iii] = ( *iter_trk )->getLErrorE()[jj][kk];
5300 m_lerrormu[iii] = ( *iter_trk )->getLErrorMu()[jj][kk];
5301 m_lerrork[iii] = ( *iter_trk )->getLErrorK()[jj][kk];
5302 m_lerrorp[iii] = ( *iter_trk )->getLErrorP()[jj][kk];
5308 m_chisq[0][0] = ( *iter_trk )->getChisq( 0, 0 );
5309 m_chisq[1][0] = ( *iter_trk )->getChisq( 0, 1 );
5310 m_chisq[2][0] = ( *iter_trk )->getChisq( 0, 2 );
5311 m_chisq[3][0] = ( *iter_trk )->getChisq( 0, 3 );
5312 m_chisq[4][0] = ( *iter_trk )->getChisq( 0, 4 );
5313 m_chisq[0][1] = ( *iter_trk )->getChisq( 1, 0 );
5314 m_chisq[1][1] = ( *iter_trk )->getChisq( 1, 1 );
5315 m_chisq[2][1] = ( *iter_trk )->getChisq( 1, 2 );
5316 m_chisq[3][1] = ( *iter_trk )->getChisq( 1, 3 );
5317 m_chisq[4][1] = ( *iter_trk )->getChisq( 1, 4 );
5319 m_ndf[0][0] = ( *iter_trk )->getNdf( 0, 0 );
5320 m_ndf[1][0] = ( *iter_trk )->getNdf( 0, 1 );
5321 m_ndf[2][0] = ( *iter_trk )->getNdf( 0, 2 );
5322 m_ndf[3][0] = ( *iter_trk )->getNdf( 0, 3 );
5323 m_ndf[4][0] = ( *iter_trk )->getNdf( 0, 4 );
5324 m_ndf[0][1] = ( *iter_trk )->getNdf( 1, 0 );
5325 m_ndf[1][1] = ( *iter_trk )->getNdf( 1, 1 );
5326 m_ndf[2][1] = ( *iter_trk )->getNdf( 1, 2 );
5327 m_ndf[3][1] = ( *iter_trk )->getNdf( 1, 3 );
5328 m_ndf[4][1] = ( *iter_trk )->getNdf( 1, 4 );
5330 m_stat[0][0] = ( *iter_trk )->getStat( 0, 0 );
5331 m_stat[1][0] = ( *iter_trk )->getStat( 0, 1 );
5332 m_stat[2][0] = ( *iter_trk )->getStat( 0, 2 );
5333 m_stat[3][0] = ( *iter_trk )->getStat( 0, 3 );
5334 m_stat[4][0] = ( *iter_trk )->getStat( 0, 4 );
5335 m_stat[0][1] = ( *iter_trk )->getStat( 1, 0 );
5336 m_stat[1][1] = ( *iter_trk )->getStat( 1, 1 );
5337 m_stat[2][1] = ( *iter_trk )->getStat( 1, 2 );
5338 m_stat[3][1] = ( *iter_trk )->getStat( 1, 3 );
5339 m_stat[4][1] = ( *iter_trk )->getStat( 1, 4 );
5341 m_fptot = sqrt( 1 + pow( m_fhelix[4], 2 ) ) / m_fhelix[2];
5342 m_fptote = sqrt( 1 + pow( m_fhelixe[4], 2 ) ) / m_fhelixe[2];
5343 m_fptotmu = sqrt( 1 + pow( m_fhelixmu[4], 2 ) ) / m_fhelixmu[2];
5344 m_fptotk = sqrt( 1 + pow( m_fhelixk[4], 2 ) ) / m_fhelixk[2];
5345 m_fptotp = sqrt( 1 + pow( m_fhelixp[4], 2 ) ) / m_fhelixp[2];
5347 m_zpt = 1 / m_zhelix[2];
5348 m_zpte = 1 / m_zhelixe[2];
5349 m_zptmu = 1 / m_zhelixmu[2];
5350 m_zptk = 1 / m_zhelixk[2];
5351 m_zptp = 1 / m_zhelixp[2];
5353 m_fpt = 1 / m_fhelix[2];
5354 m_fpte = 1 / m_fhelixe[2];
5355 m_fptmu = 1 / m_fhelixmu[2];
5356 m_fptk = 1 / m_fhelixk[2];
5357 m_fptp = 1 / m_fhelixp[2];
5359 m_lpt = 1 / m_lhelix[2];
5360 m_lpte = 1 / m_lhelixe[2];
5361 m_lptmu = 1 / m_lhelixmu[2];
5362 m_lptk = 1 / m_lhelixk[2];
5363 m_lptp = 1 / m_lhelixp[2];
5365 m_lptot = sqrt( 1 + pow( m_lhelix[4], 2 ) ) / m_lhelix[2];
5366 m_lptote = sqrt( 1 + pow( m_lhelixe[4], 2 ) ) / m_lhelixe[2];
5367 m_lptotmu = sqrt( 1 + pow( m_lhelixmu[4], 2 ) ) / m_lhelixmu[2];
5368 m_lptotk = sqrt( 1 + pow( m_lhelixk[4], 2 ) ) / m_lhelixk[2];
5369 m_lptotp = sqrt( 1 + pow( m_lhelixp[4], 2 ) ) / m_lhelixp[2];
5371 m_zptot = sqrt( 1 + pow( m_zhelix[4], 2 ) ) / m_zhelix[2];
5372 m_zptote = sqrt( 1 + pow( m_zhelixe[4], 2 ) ) / m_zhelixe[2];
5373 m_zptotmu = sqrt( 1 + pow( m_zhelixmu[4], 2 ) ) / m_zhelixmu[2];
5374 m_zptotk = sqrt( 1 + pow( m_zhelixk[4], 2 ) ) / m_zhelixk[2];
5375 m_zptotp = sqrt( 1 + pow( m_zhelixp[4], 2 ) ) / m_zhelixp[2];
5378 m_zsigp = sqrt( pow( ( m_zptot / m_zhelix[2] ), 2 ) * m_zerror[5] +
5379 pow( ( m_zhelix[4] / m_zptot ), 2 ) * pow( ( 1 / m_zhelix[2] ), 4 ) *
5381 2 * m_zhelix[4] * m_zerror[12] * pow( ( 1 / m_zhelix[2] ), 3 ) );
5382 m_zsigpe = sqrt( pow( ( m_zptote / m_zhelixe[2] ), 2 ) * m_zerrore[5] +
5383 pow( ( m_zhelixe[4] / m_zptote ), 2 ) *
5384 pow( ( 1 / m_zhelixe[2] ), 4 ) * m_zerrore[14] -
5385 2 * m_zhelixe[4] * m_zerrore[12] * pow( ( 1 / m_zhelixe[2] ), 3 ) );
5387 sqrt( pow( ( m_zptotmu / m_zhelixmu[2] ), 2 ) * m_zerrormu[5] +
5388 pow( ( m_zhelixmu[4] / m_zptotmu ), 2 ) * pow( ( 1 / m_zhelixmu[2] ), 4 ) *
5390 2 * m_zhelixmu[4] * m_zerrormu[12] * pow( ( 1 / m_zhelixmu[2] ), 3 ) );
5391 m_zsigpk = sqrt( pow( ( m_zptotk / m_zhelixk[2] ), 2 ) * m_zerrork[5] +
5392 pow( ( m_zhelixk[4] / m_zptotk ), 2 ) *
5393 pow( ( 1 / m_zhelixk[2] ), 4 ) * m_zerrork[14] -
5394 2 * m_zhelixk[4] * m_zerrork[12] * pow( ( 1 / m_zhelixk[2] ), 3 ) );
5395 m_zsigpp = sqrt( pow( ( m_zptotp / m_zhelixp[2] ), 2 ) * m_zerrorp[5] +
5396 pow( ( m_zhelixp[4] / m_zptotp ), 2 ) *
5397 pow( ( 1 / m_zhelixp[2] ), 4 ) * m_zerrorp[14] -
5398 2 * m_zhelixp[4] * m_zerrorp[12] * pow( ( 1 / m_zhelixp[2] ), 3 ) );
5401 StatusCode sc1 = m_nt1->write();
5402 if ( sc1.isFailure() ) cout <<
"Ntuple1 filling failed!" << endl;
5409 phi1 = ( *iter_trk )->getFFi0();
5410 r1 = ( *iter_trk )->getFDr();
5411 z1 = ( *iter_trk )->getFDz();
5412 kap1 = ( *iter_trk )->getFCpa();
5413 tanl1 = ( *iter_trk )->getFTanl();
5416 p1 = sqrt( 1 + tanl1 * tanl1 ) / kap1;
5417 the1 =
M_PI / 2 - atan( tanl1 );
5421 phi2 = ( *iter_trk )->getFFi0();
5422 r2 = ( *iter_trk )->getFDr();
5423 z2 = ( *iter_trk )->getFDz();
5424 kap2 = ( *iter_trk )->getFCpa();
5425 tanl2 = ( *iter_trk )->getFTanl();
5428 p2 = sqrt( 1 + tanl2 * tanl2 ) / kap1;
5429 the2 =
M_PI / 2 - atan( tanl2 );
5438 m_delthe = the1 + the2;
5441 StatusCode sc2 = m_nt2->write();
5442 if ( sc2.isFailure() ) cout <<
"Ntuple2 filling failed!" << endl;
5449 if (
debug_ == 4 ) cout <<
"Kalfitting finished " << std::endl;