381 NTuplePtr nt1(
ntupleSvc(),
"FILE104/n101");
383 if ( nt1 ) m_nt1 = nt1;
385 m_nt1=
ntupleSvc()->book(
"FILE104/n101",CLID_ColumnWiseTuple,
"KalFit");
388 status = m_nt1->addItem(
"trackid",m_trackid);
389 status = m_nt1->addItem(
"stat",5,2,m_stat);
390 status = m_nt1->addItem(
"ndf",5,2,m_ndf);
391 status = m_nt1->addItem(
"chisq",5,2,m_chisq);
392 status = m_nt1->addItem(
"length",5,m_length);
393 status = m_nt1->addItem(
"tof",5,m_tof);
394 status = m_nt1->addItem(
"nhits",5,m_nhits);
395 status = m_nt1->addItem(
"zhelix",5,m_zhelix);
396 status = m_nt1->addItem(
"zhelixe",5,m_zhelixe);
397 status = m_nt1->addItem(
"zhelixmu",5,m_zhelixmu);
398 status = m_nt1->addItem(
"zhelixk",5,m_zhelixk);
399 status = m_nt1->addItem(
"zhelixp",5,m_zhelixp);
400 status = m_nt1->addItem(
"zptot",m_zptot);
401 status = m_nt1->addItem(
"zptote",m_zptote);
402 status = m_nt1->addItem(
"zptotmu",m_zptotmu);
403 status = m_nt1->addItem(
"zptotk",m_zptotk);
404 status = m_nt1->addItem(
"zptotp",m_zptotp);
406 status = m_nt1->addItem(
"zpt",m_zpt);
407 status = m_nt1->addItem(
"zpte",m_zpte);
408 status = m_nt1->addItem(
"zptmu",m_zptmu);
409 status = m_nt1->addItem(
"zptk",m_zptk);
410 status = m_nt1->addItem(
"zptp",m_zptp);
412 status = m_nt1->addItem(
"fptot",m_fptot);
413 status = m_nt1->addItem(
"fptote",m_fptote);
414 status = m_nt1->addItem(
"fptotmu",m_fptotmu);
415 status = m_nt1->addItem(
"fptotk",m_fptotk);
416 status = m_nt1->addItem(
"fptotp",m_fptotp);
417 status = m_nt1->addItem(
"fpt",m_fpt);
418 status = m_nt1->addItem(
"fpte",m_fpte);
419 status = m_nt1->addItem(
"fptmu",m_fptmu);
420 status = m_nt1->addItem(
"fptk",m_fptk);
421 status = m_nt1->addItem(
"fptp",m_fptp);
423 status = m_nt1->addItem(
"lptot",m_lptot);
424 status = m_nt1->addItem(
"lptote",m_lptote);
425 status = m_nt1->addItem(
"lptotmu",m_lptotmu);
426 status = m_nt1->addItem(
"lptotk",m_lptotk);
427 status = m_nt1->addItem(
"lptotp",m_lptotp);
428 status = m_nt1->addItem(
"lpt",m_lpt);
429 status = m_nt1->addItem(
"lpte",m_lpte);
430 status = m_nt1->addItem(
"lptmu",m_lptmu);
431 status = m_nt1->addItem(
"lptk",m_lptk);
432 status = m_nt1->addItem(
"lptp",m_lptp);
434 status = m_nt1->addItem(
"zsigp",m_zsigp);
435 status = m_nt1->addItem(
"zsigpe",m_zsigpe);
436 status = m_nt1->addItem(
"zsigpmu",m_zsigpmu);
437 status = m_nt1->addItem(
"zsigpk",m_zsigpk);
438 status = m_nt1->addItem(
"zsigpp",m_zsigpp);
439 status = m_nt1->addItem(
"fhelix",5,m_fhelix);
440 status = m_nt1->addItem(
"fhelixe",5,m_fhelixe);
441 status = m_nt1->addItem(
"fhelixmu",5,m_fhelixmu);
442 status = m_nt1->addItem(
"fhelixk",5,m_fhelixk);
443 status = m_nt1->addItem(
"fhelixp",5,m_fhelixp);
444 status = m_nt1->addItem(
"lhelix",5,m_lhelix);
445 status = m_nt1->addItem(
"lhelixe",5,m_lhelixe);
446 status = m_nt1->addItem(
"lhelixmu",5,m_lhelixmu);
447 status = m_nt1->addItem(
"lhelixk",5,m_lhelixk);
448 status = m_nt1->addItem(
"lhelixp",5,m_lhelixp);
449 status = m_nt1->addItem(
"rGem",4,m_rGem);
450 status = m_nt1->addItem(
"chi2Gem",4,m_chi2Gem);
451 status = m_nt1->addItem(
"phiGemExp",4,m_phiGemExp);
452 status = m_nt1->addItem(
"phiGemHit",4,m_phiGemHit);
453 status = m_nt1->addItem(
"zGemExp",4,m_zGemExp);
454 status = m_nt1->addItem(
"zGemHit",4,m_zGemHit);
456 status = m_nt1->addItem(
"zerror",15,m_zerror);
457 status = m_nt1->addItem(
"zerrore",15,m_zerrore);
458 status = m_nt1->addItem(
"zerrormu",15,m_zerrormu);
459 status = m_nt1->addItem(
"zerrork",15,m_zerrork);
460 status = m_nt1->addItem(
"zerrorp",15,m_zerrorp);
461 status = m_nt1->addItem(
"ferror",15,m_ferror);
462 status = m_nt1->addItem(
"ferrore",15,m_ferrore);
463 status = m_nt1->addItem(
"ferrormu",15,m_ferrormu);
464 status = m_nt1->addItem(
"ferrork",15,m_ferrork);
465 status = m_nt1->addItem(
"ferrorp",15,m_ferrorp);
466 status = m_nt1->addItem(
"lerror",15,m_lerror);
467 status = m_nt1->addItem(
"lerrore",15,m_lerrore);
468 status = m_nt1->addItem(
"lerrormu",15,m_lerrormu);
469 status = m_nt1->addItem(
"lerrork",15,m_lerrork);
470 status = m_nt1->addItem(
"lerrorp",15,m_lerrorp);
473 status = m_nt1->addItem(
"evtid",m_evtid);
474 status = m_nt1->addItem(
"mchelix",5,m_mchelix);
475 status = m_nt1->addItem(
"mcptot",m_mcptot);
476 status = m_nt1->addItem(
"mcpid",m_mcpid);
478 if( status.isFailure() ) cout<<
"Ntuple1 add item failed!"<<endl;
484 NTuplePtr nt2(
ntupleSvc(),
"FILE104/n102");
486 if ( nt2 ) m_nt2 = nt2;
488 m_nt2=
ntupleSvc()->book(
"FILE104/n102",CLID_ColumnWiseTuple,
"KalFitComp");
490 status2 = m_nt2->addItem(
"delx",m_delx);
491 status2 = m_nt2->addItem(
"dely",m_dely);
492 status2 = m_nt2->addItem(
"delz",m_delz);
493 status2 = m_nt2->addItem(
"delthe",m_delthe);
494 status2 = m_nt2->addItem(
"delphi",m_delphi);
495 status2 = m_nt2->addItem(
"delp",m_delp);
496 status2 = m_nt2->addItem(
"delpx",m_delpx);
497 status2 = m_nt2->addItem(
"delpy",m_delpy);
498 status2 = m_nt2->addItem(
"delpz",m_delpz);
500 if( status2.isFailure() ) cout<<
"Ntuple2 add item failed!"<<endl;
506 NTuplePtr nt3(
ntupleSvc(),
"FILE104/n103");
508 if ( nt3 ) m_nt3 = nt3;
510 m_nt3=
ntupleSvc()->book(
"FILE104/n103",CLID_ColumnWiseTuple,
"PatRec");
512 status3 = m_nt3->addItem(
"trkhelix",5,m_trkhelix);
513 status3 = m_nt3->addItem(
"trkptot",m_trkptot);
515 status3 = m_nt3->addItem(
"trkerror",15,m_trkerror);
516 status3 = m_nt3->addItem(
"trksigp",m_trksigp);
518 status3 = m_nt3->addItem(
"trkndf",m_trkndf);
519 status3 = m_nt3->addItem(
"trkchisq",m_trkchisq);
520 if( status3.isFailure() ) cout<<
"Ntuple3 add item failed!"<<endl;
526 NTuplePtr nt4(
ntupleSvc(),
"FILE104/n104");
528 if ( nt4 ) m_nt4 = nt4;
530 m_nt4=
ntupleSvc()->book(
"FILE104/n104",CLID_ColumnWiseTuple,
"PatRecComp");
532 status4 = m_nt4->addItem(
"trkdelx",m_trkdelx);
533 status4 = m_nt4->addItem(
"trkdely",m_trkdely);
534 status4 = m_nt4->addItem(
"trkdelz",m_trkdelz);
535 status4 = m_nt4->addItem(
"trkdelthe",m_trkdelthe);
536 status4 = m_nt4->addItem(
"trkdelphi",m_trkdelphi);
537 status4 = m_nt4->addItem(
"trkdelp",m_trkdelp);
538 if( status4.isFailure() ) cout<<
"Ntuple4 add item failed!"<<endl;
543 NTuplePtr nt5(
ntupleSvc(),
"FILE104/n105");
545 if ( nt5 ) m_nt5 = nt5;
547 m_nt5=
ntupleSvc()->book(
"FILE104/n105",CLID_ColumnWiseTuple,
"KalFitdChisq");
549 status5 = m_nt5->addItem(
"dchi2",m_dchi2);
550 status5 = m_nt5->addItem(
"masshyp",m_masshyp);
551 status5 = m_nt5->addItem(
"residual_estim",m_residest);
552 status5 = m_nt5->addItem(
"residual",m_residnew);
553 status5 = m_nt5->addItem(
"layer",m_layer);
554 status5 = m_nt5->addItem(
"kaldr",m_anal_dr);
555 status5 = m_nt5->addItem(
"kalphi0",m_anal_phi0);
556 status5 = m_nt5->addItem(
"kalkappa",m_anal_kappa);
557 status5 = m_nt5->addItem(
"kaldz",m_anal_dz);
558 status5 = m_nt5->addItem(
"kaltanl",m_anal_tanl);
559 status5 = m_nt5->addItem(
"dr_ea",m_anal_ea_dr);
560 status5 = m_nt5->addItem(
"phi0_ea",m_anal_ea_phi0);
561 status5 = m_nt5->addItem(
"kappa_ea",m_anal_ea_kappa);
562 status5 = m_nt5->addItem(
"dz_ea",m_anal_ea_dz);
563 status5 = m_nt5->addItem(
"tanl_ea",m_anal_ea_tanl);
564 if( status5.isFailure() ) cout<<
"Ntuple5 add item failed!"<<endl;
571 NTuplePtr nt7(
ntupleSvc(),
"FILE104/n120");
575 m_nt7 =
ntupleSvc()-> book(
"FILE104/n120",CLID_ColumnWiseTuple,
"kalfit_failure");
577 status7 = m_nt7 ->addItem(
"run_kal",m_run_kal);
578 status7 = m_nt7 ->addItem(
"event_kal",m_event_kal);
579 status7 = m_nt7 ->addItem(
"trkid_kal",m_trkid_kal);
580 status7 = m_nt7 ->addItem(
"dropedHits_kal_e",m_dropedHits_kal_e);
581 status7 = m_nt7 ->addItem(
"kappa2_kal_e",m_kappa2_kal_e);
582 status7 = m_nt7 ->addItem(
"trackNhits_kal_e",m_trackNhits_kal_e);
583 status7 = m_nt7 ->addItem(
"trackNster_kal_e",m_trackNster_kal_e);
584 status7 = m_nt7 ->addItem(
"trackNaxis_kal_e",m_trackNaxis_kal_e);
585 status7 = m_nt7 ->addItem(
"chi2_kal_e",m_chi2_kal_e);
586 status7 = m_nt7 ->addItem(
"Ea00_kal_e",m_Ea00_kal_e);
587 status7 = m_nt7 ->addItem(
"Ea11_kal_e",m_Ea11_kal_e);
588 status7 = m_nt7 ->addItem(
"Ea22_kal_e",m_Ea22_kal_e);
589 status7 = m_nt7 ->addItem(
"Ea33_kal_e",m_Ea33_kal_e);
590 status7 = m_nt7 ->addItem(
"Ea44_kal_e",m_Ea44_kal_e);
591 status7 = m_nt7 ->addItem(
"dropedHits_kal_mu",m_dropedHits_kal_mu);
592 status7 = m_nt7 ->addItem(
"kappa2_kal_mu",m_kappa2_kal_mu);
593 status7 = m_nt7 ->addItem(
"trackNhits_kal_mu",m_trackNhits_kal_mu);
594 status7 = m_nt7 ->addItem(
"trackNster_kal_mu",m_trackNster_kal_mu);
595 status7 = m_nt7 ->addItem(
"trackNaxis_kal_mu",m_trackNaxis_kal_mu);
596 status7 = m_nt7 ->addItem(
"chi2_kal_mu",m_chi2_kal_mu);
597 status7 = m_nt7 ->addItem(
"Ea00_kal_mu",m_Ea00_kal_mu);
598 status7 = m_nt7 ->addItem(
"Ea11_kal_mu",m_Ea11_kal_mu);
599 status7 = m_nt7 ->addItem(
"Ea22_kal_mu",m_Ea22_kal_mu);
600 status7 = m_nt7 ->addItem(
"Ea33_kal_mu",m_Ea33_kal_mu);
601 status7 = m_nt7 ->addItem(
"Ea44_kal_mu",m_Ea44_kal_mu);
602 status7 = m_nt7 ->addItem(
"iqual_front_kal_mu",m_iqual_front_kal_mu);
603 status7 = m_nt7 ->addItem(
"dropedHits_kal_pi",m_dropedHits_kal_pi);
604 status7 = m_nt7 ->addItem(
"kappa2_kal_pi",m_kappa2_kal_pi);
605 status7 = m_nt7 ->addItem(
"trackNhits_kal_pi",m_trackNhits_kal_pi);
606 status7 = m_nt7 ->addItem(
"trackNster_kal_pi",m_trackNster_kal_pi);
607 status7 = m_nt7 ->addItem(
"trackNaxis_kal_pi",m_trackNaxis_kal_pi);
608 status7 = m_nt7 ->addItem(
"chi2_kal_pi",m_chi2_kal_pi);
609 status7 = m_nt7 ->addItem(
"Ea00_kal_pi",m_Ea00_kal_pi);
610 status7 = m_nt7 ->addItem(
"Ea11_kal_pi",m_Ea11_kal_pi);
611 status7 = m_nt7 ->addItem(
"Ea22_kal_pi",m_Ea22_kal_pi);
612 status7 = m_nt7 ->addItem(
"Ea33_kal_pi",m_Ea33_kal_pi);
613 status7 = m_nt7 ->addItem(
"Ea44_kal_pi",m_Ea44_kal_pi);
614 status7 = m_nt7 ->addItem(
"dropedHits_kal_k",m_dropedHits_kal_k);
615 status7 = m_nt7 ->addItem(
"kappa2_kal_k",m_kappa2_kal_k);
616 status7 = m_nt7 ->addItem(
"trackNhits_kal_k",m_trackNhits_kal_k);
617 status7 = m_nt7 ->addItem(
"trackNster_kal_k",m_trackNster_kal_k);
618 status7 = m_nt7 ->addItem(
"trackNaxis_kal_k",m_trackNaxis_kal_k);
619 status7 = m_nt7 ->addItem(
"chi2_kal_k",m_chi2_kal_k);
620 status7 = m_nt7 ->addItem(
"Ea00_kal_k",m_Ea00_kal_k);
621 status7 = m_nt7 ->addItem(
"Ea11_kal_k",m_Ea11_kal_k);
622 status7 = m_nt7 ->addItem(
"Ea22_kal_k",m_Ea22_kal_k);
623 status7 = m_nt7 ->addItem(
"Ea33_kal_k",m_Ea33_kal_k);
624 status7 = m_nt7 ->addItem(
"Ea44_kal_k",m_Ea44_kal_k);
625 status7 = m_nt7 ->addItem(
"dropedHits_kal_p",m_dropedHits_kal_p);
626 status7 = m_nt7 ->addItem(
"kappa2_kal_p",m_kappa2_kal_p);
627 status7 = m_nt7 ->addItem(
"trackNhits_kal_p",m_trackNhits_kal_p);
628 status7 = m_nt7 ->addItem(
"trackNster_kal_p",m_trackNster_kal_p);
629 status7 = m_nt7 ->addItem(
"trackNaxis_kal_p",m_trackNaxis_kal_p);
630 status7 = m_nt7 ->addItem(
"chi2_kal_p",m_chi2_kal_p);
631 status7 = m_nt7 ->addItem(
"Ea00_kal_p",m_Ea00_kal_p);
632 status7 = m_nt7 ->addItem(
"Ea11_kal_p",m_Ea11_kal_p);
633 status7 = m_nt7 ->addItem(
"Ea22_kal_p",m_Ea22_kal_p);
634 status7 = m_nt7 ->addItem(
"Ea33_kal_p",m_Ea33_kal_p);
635 status7 = m_nt7 ->addItem(
"Ea44_kal_p",m_Ea44_kal_p);
637 status7 = m_nt7 ->addItem(
"hit_number",m_hit_no,0, 1000);
638 status7 = m_nt7 ->addItem(
"nCluster",m_nCluster,0, 1000);
639 status7 = m_nt7->addIndexedItem(
"dchi2_hit_e",m_hit_no,m_dchi2_hit_e);
640 status7 = m_nt7->addIndexedItem(
"residual_estim_hit_e",m_hit_no,m_residest_hit_e);
641 status7 = m_nt7->addIndexedItem(
"residual_hit_e",m_hit_no,m_residnew_hit_e);
642 status7 = m_nt7->addIndexedItem(
"layer_hit_e",m_hit_no,m_layer_hit_e);
643 status7 = m_nt7->addIndexedItem(
"kaldr_hit_e",m_hit_no,m_anal_dr_hit_e);
644 status7 = m_nt7->addIndexedItem(
"kalphi0_hit_e",m_hit_no,m_anal_phi0_hit_e);
645 status7 = m_nt7->addIndexedItem(
"kalkappa_hit_e",m_hit_no,m_anal_kappa_hit_e);
646 status7 = m_nt7->addIndexedItem(
"kaldz_hit_e",m_hit_no,m_anal_dz_hit_e);
647 status7 = m_nt7->addIndexedItem(
"kaltanl_hit_e",m_hit_no,m_anal_tanl_hit_e);
648 status7 = m_nt7->addIndexedItem(
"dr_ea_hit_e",m_hit_no,m_anal_ea_dr_hit_e);
649 status7 = m_nt7->addIndexedItem(
"phi0_ea_hit_e",m_hit_no,m_anal_ea_phi0_hit_e);
650 status7 = m_nt7->addIndexedItem(
"kappa_ea_hit_e",m_hit_no,m_anal_ea_kappa_hit_e);
651 status7 = m_nt7->addIndexedItem(
"dz_ea_hit_e",m_hit_no,m_anal_ea_dz_hit_e);
652 status7 = m_nt7->addIndexedItem(
"tanl_ea_hit_e",m_hit_no,m_anal_ea_tanl_hit_e);
653 status7 = m_nt7->addIndexedItem(
"dchi2_hit_mu",m_hit_no,m_dchi2_hit_mu);
654 status7 = m_nt7->addIndexedItem(
"residual_estim_hit_mu",m_hit_no,m_residest_hit_mu);
655 status7 = m_nt7->addIndexedItem(
"residual_hit_mu",m_hit_no,m_residnew_hit_mu);
656 status7 = m_nt7->addIndexedItem(
"layer_hit_mu",m_hit_no,m_layer_hit_mu);
657 status7 = m_nt7->addIndexedItem(
"kaldr_hit_mu",m_hit_no,m_anal_dr_hit_mu);
658 status7 = m_nt7->addIndexedItem(
"kalphi0_hit_mu",m_hit_no,m_anal_phi0_hit_mu);
659 status7 = m_nt7->addIndexedItem(
"kalkappa_hit_mu",m_hit_no,m_anal_kappa_hit_mu);
660 status7 = m_nt7->addIndexedItem(
"kaldz_hit_mu",m_hit_no,m_anal_dz_hit_mu);
661 status7 = m_nt7->addIndexedItem(
"kaltanl_hit_mu",m_hit_no,m_anal_tanl_hit_mu);
662 status7 = m_nt7->addIndexedItem(
"dr_ea_hit_mu",m_hit_no,m_anal_ea_dr_hit_mu);
663 status7 = m_nt7->addIndexedItem(
"phi0_ea_hit_mu",m_hit_no,m_anal_ea_phi0_hit_mu);
664 status7 = m_nt7->addIndexedItem(
"kappa_ea_hit_mu",m_hit_no,m_anal_ea_kappa_hit_mu);
665 status7 = m_nt7->addIndexedItem(
"dz_ea_hit_mu",m_hit_no,m_anal_ea_dz_hit_mu);
666 status7 = m_nt7->addIndexedItem(
"tanl_ea_hit_mu",m_hit_no,m_anal_ea_tanl_hit_mu);
667 status7 = m_nt7->addIndexedItem(
"dchi2_hit_pi",m_hit_no,m_dchi2_hit_pi);
668 status7 = m_nt7->addIndexedItem(
"residual_estim_hit_pi",m_hit_no,m_residest_hit_pi);
669 status7 = m_nt7->addIndexedItem(
"residual_hit_pi",m_hit_no,m_residnew_hit_pi);
670 status7 = m_nt7->addIndexedItem(
"layer_hit_pi",m_hit_no,m_layer_hit_pi);
671 status7 = m_nt7->addIndexedItem(
"kaldr_hit_pi",m_hit_no,m_anal_dr_hit_pi);
672 status7 = m_nt7->addIndexedItem(
"kalphi0_hit_pi",m_hit_no,m_anal_phi0_hit_pi);
673 status7 = m_nt7->addIndexedItem(
"kalkappa_hit_pi",m_hit_no,m_anal_kappa_hit_pi);
674 status7 = m_nt7->addIndexedItem(
"kaldz_hit_pi",m_hit_no,m_anal_dz_hit_pi);
675 status7 = m_nt7->addIndexedItem(
"kaltanl_hit_pi",m_hit_no,m_anal_tanl_hit_pi);
676 status7 = m_nt7->addIndexedItem(
"dr_ea_hit_pi",m_hit_no,m_anal_ea_dr_hit_pi);
677 status7 = m_nt7->addIndexedItem(
"phi0_ea_hit_pi",m_hit_no,m_anal_ea_phi0_hit_pi);
678 status7 = m_nt7->addIndexedItem(
"kappa_ea_hit_pi",m_hit_no,m_anal_ea_kappa_hit_pi);
679 status7 = m_nt7->addIndexedItem(
"dz_ea_hit_pi",m_hit_no,m_anal_ea_dz_hit_pi);
680 status7 = m_nt7->addIndexedItem(
"tanl_ea_hit_pi",m_hit_no,m_anal_ea_tanl_hit_pi);
681 status7 = m_nt7->addIndexedItem(
"dchi2_hit_k",m_hit_no,m_dchi2_hit_k);
682 status7 = m_nt7->addIndexedItem(
"residual_estim_hit_k",m_hit_no,m_residest_hit_k);
683 status7 = m_nt7->addIndexedItem(
"residual_hit_k",m_hit_no,m_residnew_hit_k);
684 status7 = m_nt7->addIndexedItem(
"layer_hit_k",m_hit_no,m_layer_hit_k);
685 status7 = m_nt7->addIndexedItem(
"kaldr_hit_k",m_hit_no,m_anal_dr_hit_k);
686 status7 = m_nt7->addIndexedItem(
"kalphi0_hit_k",m_hit_no,m_anal_phi0_hit_k);
687 status7 = m_nt7->addIndexedItem(
"kalkappa_hit_k",m_hit_no,m_anal_kappa_hit_k);
688 status7 = m_nt7->addIndexedItem(
"kaldz_hit_k",m_hit_no,m_anal_dz_hit_k);
689 status7 = m_nt7->addIndexedItem(
"kaltanl_hit_k",m_hit_no,m_anal_tanl_hit_k);
690 status7 = m_nt7->addIndexedItem(
"dr_ea_hit_k",m_hit_no,m_anal_ea_dr_hit_k);
691 status7 = m_nt7->addIndexedItem(
"phi0_ea_hit_k",m_hit_no,m_anal_ea_phi0_hit_k);
692 status7 = m_nt7->addIndexedItem(
"kappa_ea_hit_k",m_hit_no,m_anal_ea_kappa_hit_k);
693 status7 = m_nt7->addIndexedItem(
"dz_ea_hit_k",m_hit_no,m_anal_ea_dz_hit_k);
694 status7 = m_nt7->addIndexedItem(
"tanl_ea_hit_k",m_hit_no,m_anal_ea_tanl_hit_k);
695 status7 = m_nt7->addIndexedItem(
"dchi2_hit_p",m_hit_no,m_dchi2_hit_p);
696 status7 = m_nt7->addIndexedItem(
"residual_estim_hit_p",m_hit_no,m_residest_hit_p);
697 status7 = m_nt7->addIndexedItem(
"residual_hit_p",m_hit_no,m_residnew_hit_p);
698 status7 = m_nt7->addIndexedItem(
"layer_hit_p",m_hit_no,m_layer_hit_p);
699 status7 = m_nt7->addIndexedItem(
"kaldr_hit_p",m_hit_no,m_anal_dr_hit_p);
700 status7 = m_nt7->addIndexedItem(
"kalphi0_hit_p",m_hit_no,m_anal_phi0_hit_p);
701 status7 = m_nt7->addIndexedItem(
"kalkappa_hit_p",m_hit_no,m_anal_kappa_hit_p);
702 status7 = m_nt7->addIndexedItem(
"kaldz_hit_p",m_hit_no,m_anal_dz_hit_p);
703 status7 = m_nt7->addIndexedItem(
"kaltanl_hit_p",m_hit_no,m_anal_tanl_hit_p);
704 status7 = m_nt7->addIndexedItem(
"dr_ea_hit_p",m_hit_no,m_anal_ea_dr_hit_p);
705 status7 = m_nt7->addIndexedItem(
"phi0_ea_hit_p",m_hit_no,m_anal_ea_phi0_hit_p);
706 status7 = m_nt7->addIndexedItem(
"kappa_ea_hit_p",m_hit_no,m_anal_ea_kappa_hit_p);
707 status7 = m_nt7->addIndexedItem(
"dz_ea_hit_p",m_hit_no,m_anal_ea_dz_hit_p);
708 status7 = m_nt7->addIndexedItem(
"tanl_ea_hit_p",m_hit_no,m_anal_ea_tanl_hit_p);
710 if( status7.isFailure() ) cout<<
"Ntuple7 add item failed!"<<endl;
717 NTuplePtr nt6(
ntupleSvc(),
"FILE104/n106");
719 if ( nt6 ) m_nt6 = nt6;
721 m_nt6=
ntupleSvc()->book(
"FILE104/n106",CLID_ColumnWiseTuple,
"kal seg");
723 status6 = m_nt6->addItem(
"docaInc",m_docaInc);
724 status6 = m_nt6->addItem(
"docaExc",m_docaExc);
725 status6 = m_nt6->addItem(
"tdr",m_tdrift);
726 status6 = m_nt6->addItem(
"layerid", m_layerid);
727 status6 = m_nt6->addItem(
"event", m_eventNo);
728 status6 = m_nt6->addItem(
"residualInc",m_residualInc);
729 status6 = m_nt6->addItem(
"residualExc",m_residualExc);
730 status6 = m_nt6->addItem(
"lr",m_lr);
731 status6 = m_nt6->addItem(
"dd",m_dd);
732 status6 = m_nt6->addItem(
"y",m_yposition);
734 if( status6.isFailure() ) cout<<
"Ntuple6 add item failed!"<<endl;
739 NTuplePtr nt10(
ntupleSvc(),
"FILE104/n110");
741 if ( nt10 ) m_nt10 = nt10;
743 m_nt10=
ntupleSvc()->book(
"FILE104/n110",CLID_ColumnWiseTuple,
"test");
745 status10 = m_nt10->addItem(
"evt",m_evt3);
746 status10 = m_nt10->addItem(
"qua",m_qua);
747 status10 = m_nt10->addItem(
"nhit",m_nGemHits,0,4000000);
748 status10 = m_nt10->addIndexedItem(
"meas_r",m_nGemHits,m_meas_r);
749 status10 = m_nt10->addIndexedItem(
"meas_phi",m_nGemHits,m_meas_phi);
750 status10 = m_nt10->addIndexedItem(
"meas_z",m_nGemHits,m_meas_z);
751 status10 = m_nt10->addIndexedItem(
"esti1_r",m_nGemHits,m_esti1_r);
752 status10 = m_nt10->addIndexedItem(
"esti1_phi",m_nGemHits,m_esti1_phi);
753 status10 = m_nt10->addIndexedItem(
"esti1_z",m_nGemHits,m_esti1_z);
754 status10 = m_nt10->addIndexedItem(
"esti2_r",m_nGemHits,m_esti2_r);
755 status10 = m_nt10->addIndexedItem(
"esti2_phi",m_nGemHits,m_esti2_phi);
756 status10 = m_nt10->addIndexedItem(
"esti2_z",m_nGemHits,m_esti2_z);
757 status10 = m_nt10->addIndexedItem(
"diff1_phi",m_nGemHits,m_diff1_phi);
758 status10 = m_nt10->addIndexedItem(
"diff1_z",m_nGemHits,m_diff1_z);
759 status10 = m_nt10->addIndexedItem(
"diff2_phi",m_nGemHits,m_diff2_phi);
760 status10 = m_nt10->addIndexedItem(
"diff2_z",m_nGemHits,m_diff2_z);
761 status10 = m_nt10->addIndexedItem(
"layer",m_nGemHits,m_GemLayer);
762 status10 = m_nt10->addIndexedItem(
"mass",m_nGemHits,m_mass);
763 status10 = m_nt10->addIndexedItem(
"dchi2",m_nGemHits,m_Gchi2);
764 status10 = m_nt10->addIndexedItem(
"meas_phierr",m_nGemHits,m_meas_phierr);
765 status10 = m_nt10->addIndexedItem(
"meas_zerr",m_nGemHits,m_meas_zerr);
766 status10 = m_nt10->addIndexedItem(
"esti_phierr",m_nGemHits,m_esti_phierr);
767 status10 = m_nt10->addIndexedItem(
"esti_zerr",m_nGemHits,m_esti_zerr);
768 if( status10.isFailure() ) cout<<
"Ntuple10 add item failed!"<<endl;
774 NTuplePtr nt11(
ntupleSvc(),
"FILE104/n111");
776 if ( nt11 ) m_nt11 = nt11;
778 m_nt11=
ntupleSvc()->book(
"FILE104/n111",CLID_ColumnWiseTuple,
"truth");
780 status11 = m_nt11->addItem(
"evt",m_evt4);
781 status11 = m_nt11->addItem(
"ntruth",m_ntruth,0,400000);
782 status11 = m_nt11->addIndexedItem(
"phi",m_ntruth,m_dtphi);
783 status11 = m_nt11->addIndexedItem(
"z",m_ntruth,m_dtv);
784 status11 = m_nt11->addIndexedItem(
"postphi",m_ntruth,m_dtpostphi);
785 status11 = m_nt11->addIndexedItem(
"postz",m_ntruth,m_dtpostz);
786 status11 = m_nt11->addIndexedItem(
"layer",m_ntruth,m_tlayer);
787 if( status11.isFailure() ) cout<<
"Ntuple11 add item failed!"<<endl;
793 NTuplePtr nt12(
ntupleSvc(),
"FILE104/n112");
795 if ( nt12 ) m_nt12 = nt12;
797 m_nt12=
ntupleSvc()->book(
"FILE104/n112",CLID_ColumnWiseTuple,
"PatRecComp");
799 status12 = m_nt12->addItem(
"evt",m_evt);
800 status12 = m_nt12->addItem(
"track",m_track);
801 status12 = m_nt12->addItem(
"diff_dr",m_diff_dr);
802 status12 = m_nt12->addItem(
"diff_phi0",m_diff_phi0);
803 status12 = m_nt12->addItem(
"diff_kappa",m_diff_kappa);
804 status12 = m_nt12->addItem(
"diff_dz",m_diff_dz);
805 status12 = m_nt12->addItem(
"diff_tanl",m_diff_tanl);
806 status12 = m_nt12->addItem(
"diff_p",m_diff_p);
807 if( status12.isFailure() ) cout<<
"Ntuple12 add item failed!"<<endl;
920 MsgStream log(
msgSvc(), name());
921 log << MSG::INFO <<
"in execute()" << endreq;
922 if(
testOutput) std::cout<<
"begin to deal with EVENT ..."<<(++
eventno)<<std::endl;
968 StatusCode sc = service (
"MagneticFieldSvc",IMFSvc);
969 if(sc!=StatusCode::SUCCESS) {
970 log << MSG::ERROR <<
"Unable to open Magnetic field service"<<endreq;
979 std::cout<<
" execute, referField from MagneticFieldSvc: "<< (IMFSvc->
getReferField())*10000 <<std::endl;
986 IDataProviderSvc* evtSvc = NULL;
987 Gaudi::svcLocator()->service(
"EventDataSvc", evtSvc);
989 log << MSG::INFO <<
"makeTds:event Svc has been found" << endreq;
991 log << MSG::FATAL <<
"makeTds:Could not find eventSvc" << endreq;
992 return StatusCode::SUCCESS;
996 IDataManagerSvc *dataManSvc;
997 dataManSvc=
dynamic_cast<IDataManagerSvc*
>(evtSvc);
998 DataObject *aKalTrackCol;
999 evtSvc->findObject(
"/Event/Recon/RecMdcKalTrackCol",aKalTrackCol);
1000 if(aKalTrackCol != NULL) {
1001 dataManSvc->clearSubTree(
"/Event/Recon/RecMdcKalTrackCol");
1002 evtSvc->unregisterObject(
"/Event/Recon/RecMdcKalTrackCol");
1005 kalsc = evtSvc->registerObject(
"/Event/Recon/RecMdcKalTrackCol", kalcol);
1006 if( kalsc.isFailure() ) {
1007 log << MSG::FATAL <<
"Could not register RecMdcKalTrack" << endreq;
1008 return StatusCode::SUCCESS;
1010 log << MSG::INFO <<
"RecMdcKalTrackCol registered successfully!" <<endreq;
1012 DataObject *aKalHelixSegCol;
1013 evtSvc->findObject(
"/Event/Recon/RecMdcKalHelixSegCol", aKalHelixSegCol);
1014 if(aKalHelixSegCol != NULL){
1015 dataManSvc->clearSubTree(
"/Event/Recon/RecMdcKalHelixSegCol");
1016 evtSvc->unregisterObject(
"/Event/Recon/RecMdcKalHelixSegCol");
1019 kalsc = evtSvc->registerObject(
"/Event/Recon/RecMdcKalHelixSegCol", helixsegcol);
1020 if( kalsc.isFailure()){
1021 log<< MSG::FATAL <<
"Could not register RecMdcKalHelixSeg" <<endreq;
1022 return StatusCode::SUCCESS;
1024 log << MSG::INFO <<
"RecMdcKalHelixSegCol register successfully!" <<endreq;
1034 ISvcLocator* svcLocator = Gaudi::svcLocator();
1036 StatusCode Cgem_sc=svcLocator->service(
"CgemGeomSvc", ISvc);
1038 if (!Cgem_sc.isSuccess()) log<< MSG::INFO <<
"KalFitAlg::execute(): Could not open CGEM geometry file" << endreq;
1044 std::cout<<
"ERROR OCCUR when dynamic_cast in KalFitAlg::execute ...!!"<<std::endl;
1047 SmartDataPtr<Event::EventHeader> eventHeader(eventSvc(),
"/Event/EventHeader");
1049 log << MSG::WARNING <<
"Could not find Event Header" << endreq;
1050 return StatusCode::FAILURE;
1052 int eventNo = eventHeader->eventNumber();
1054 int runNo = eventHeader->runNumber();
1058 if(
testOutput) cout<<endl<<
"$$$$$$$$$$$ run="<<
runNo<<
", evt="<<
eventNo<<
" $$$$$$$$$$$$$$$$$"<<endl<<endl;
1061 SmartDataPtr<RecEsTimeCol> estimeCol(eventSvc(),
"/Event/Recon/RecEsTimeCol");
1062 if (estimeCol && estimeCol->size()) {
1063 RecEsTimeCol::iterator iter_evt = estimeCol->begin();
1064 t0 = (*iter_evt)->getTest();
1067 log << MSG::WARNING <<
"Could not find EvTimeCol" << endreq;
1068 return StatusCode::SUCCESS;
1073 std::cout<<
"in KalFitAlg , we get the event start time = "<<t0<<std::endl;
1077 SmartDataPtr<MdcDigiCol> mdcDigiCol(evtSvc,
"/Event/Digi/MdcDigiCol");
1078 if (sc!=StatusCode::SUCCESS) {
1079 log << MSG::FATAL <<
"Could not find MdcDigiCol!" << endreq;
1080 return StatusCode::SUCCESS;
1091 m_evtid = eventHeader->eventNumber();
1094 SmartDataPtr<McParticleCol> mcPartCol(eventSvc(),
"/Event/MC/McParticleCol");
1096 log << MSG::WARNING <<
"Could not find McParticle" << endreq;
1101 McParticleCol::iterator i_mcTrk = mcPartCol->begin();
1102 for (;i_mcTrk != mcPartCol->end(); i_mcTrk++) {
1103 if(!(*i_mcTrk)->primaryParticle())
continue;
1104 const HepLorentzVector& mom((*i_mcTrk)->initialFourMomentum());
1105 const HepLorentzVector& pos = (*i_mcTrk)->initialPosition();
1106 log << MSG::DEBUG <<
"MCINFO:particleId=" << (*i_mcTrk)->particleProperty()
1107 <<
" theta=" << mom.theta() <<
" phi="<< mom.phi()
1108 <<
" px="<< mom.px() <<
" py="<< mom.py() <<
" pz="<< mom.pz()
1110 double charge = 0.0;
1111 int pid = (*i_mcTrk)->particleProperty();
1112 int mother_id = ((*i_mcTrk)->mother()).particleProperty();
1114 charge = m_particleTable->particle( pid )->charge();
1115 }
else if ( pid <0 ) {
1116 charge = m_particleTable->particle( -pid )->charge();
1119 log << MSG::WARNING <<
"wrong particle id, please check data" <<endreq;
1122 Hep3Vector mom2(mom.px(),mom.py(),mom.pz());
1125 log << MSG::DEBUG <<
"charge of the track " << charge << endreq;
1126 if(
debug_ == 4) cout<<
"helix: "<<mchelix.
a()<<endl;
1128 for(
int j =0; j<5; j++) {
1129 m_mchelix[j] = mchelix.
a()[j];
1132 m_mcptot = sqrt(1+pow(m_mchelix[4],2))/m_mchelix[2];
1135 cout<<
"MC pid, mother_id = "<<pid<<
", "<<mother_id<<endl;
1136 cout<<
" p4 = "<<(*i_mcTrk)->initialFourMomentum()<<endl;
1137 cout<<
" start from "<<(*i_mcTrk)->initialPosition()<<endl;
1138 cout<<
" end at "<<(*i_mcTrk)->finalPosition()<<endl;
1139 cout<<
" Helix: "<<mchelix.
a()<<endl;
1140 cout<<
"mc ptot, theta, phi, R = "<<m_mcptot<<
", "<<mom2.theta()/acos(-1.)*180<<
", "<<mom2.phi()/acos(-1.)*180<<
", "<<mchelix.
radius()<<endl;
1149 SmartDataPtr<RecMdcTrackCol> newtrkCol(eventSvc(),
"/Event/Recon/RecMdcTrackCol");
1151 log << MSG::FATAL <<
"Could not find RecMdcTrackCol" << endreq;
1152 return( StatusCode::SUCCESS);
1154 log << MSG::INFO <<
"Begin to make MdcRecTrkCol and MdcRecWirhitCol"<<endreq;
1160 mtrkadd_mgr->clear();
1164 double trkx1,trkx2,trky1,trky2,trkz1,trkz2,trkthe1,trkthe2,trkphi1,trkphi2,trkp1,trkp2,trkr1,trkr2,trkkap1,trkkap2,trktanl1,trktanl2;
1165 Hep3Vector csmp3[2];
1168 RecMdcTrackCol::iterator iter_trk = newtrkCol->begin();
1169 for(
int kj = 1; iter_trk != newtrkCol->end(); iter_trk++,kj++) {
1171 csmp3[kj-1]=(*iter_trk)->p3();
1172 csmphi[kj-1] = (*iter_trk)->phi();
1176 for(
int j = 0, ij = 0; j<5; j++) {
1177 m_trkhelix[j] = (*iter_trk)->helix()[j];
1179 for(
int k=0; k<=j; k++,ij++) {
1180 m_trkerror[ij] = (*iter_trk)->err()[j][k];
1184 m_trkptot = sqrt(1+pow(m_trkhelix[4],2))/m_trkhelix[2];
1186 m_trksigp = sqrt(pow((m_trkptot/m_trkhelix[2]),2)*m_trkerror[5]+
1187 pow((m_trkhelix[4]/m_trkptot),2)*pow((1/m_trkhelix[2]),4)*m_trkerror[14]-
1188 2*m_trkhelix[4]*m_trkerror[12]*pow((1/m_trkhelix[2]),3));
1190 m_trkndf = (*iter_trk)->ndof();
1191 m_trkchisq = (*iter_trk)->chi2();
1193 if (
debug_ == 4) cout<<
"Ea from RecMdcTrackCol..." <<(*iter_trk)->err()<<endl;
1195 StatusCode sc3 = m_nt3->write();
1196 if( sc3.isFailure() ) cout<<
"Ntuple3 filling failed!"<<endl;
1226 log << MSG::DEBUG <<
"retrieved MDC tracks:"
1227 <<
" Nhits " <<(*iter_trk)->getNhits()
1228 <<
" Nster " <<(*iter_trk)->nster() <<endreq;
1230 HitRefVec gothits = (*iter_trk)->getVecHits();
1234 rectrk->
id = (*iter_trk)->trackId();
1235 rectrk->
chiSq = (*iter_trk)->chi2();
1236 rectrk->
ndf = (*iter_trk)->ndof();
1237 rectrk->
fiTerm = (*iter_trk)->getFiTerm();
1238 rectrk->
nhits = (*iter_trk)->getNhits();
1239 rectrk->
nster = (*iter_trk)->nster();
1240 rectrk->
nclus = (*iter_trk)->getNcluster();
1241 rectrk->
stat = (*iter_trk)->stat();
1242 status_temp = (*iter_trk)->stat();
1244 trkadd->
id = (*iter_trk)->trackId();
1248 trkadd->
body = rectrk;
1249 rectrk->
add = trkadd;
1254 for (
int i=0; i<5; i++) {
1255 rectrk->
helix[i] = (*iter_trk)->helix()[i];
1256 if( i<3 ) rectrk->
pivot[i] = (*iter_trk)->getPivot()[i];
1257 for(
int j = 1; j<i+2;j++) {
1258 rectrk->
error[i*(i+1)/2+j-1] = (*iter_trk)->err()(i+1,j);
1261 std::sort(gothits.begin(), gothits.end(), order_rechits);
1262 HitRefVec::iterator it_gothit = gothits.begin();
1263 for( ; it_gothit != gothits.end(); it_gothit++) {
1265 if( (*it_gothit)->getStat() != 1 ) {
1267 log<<MSG::WARNING<<
"this hit is not used in helix fitting!"<<endreq;
1272 log << MSG::DEBUG <<
"retrieved hits in MDC tracks:"
1273 <<
" hits DDL " <<(*it_gothit)->getDriftDistLeft()
1274 <<
" hits DDR " <<(*it_gothit)->getDriftDistRight()
1275 <<
" error DDL " <<(*it_gothit)->getErrDriftDistLeft()
1276 <<
" error DDR " <<(*it_gothit)->getErrDriftDistRight()
1277 <<
" id of hit "<<(*it_gothit)->getId()
1278 <<
" track id of hit "<<(*it_gothit)->getTrkId()
1279 <<
" hits ADC " <<(*it_gothit)->getAdc() << endreq;
1282 whit->
id = (*it_gothit)->getId();
1283 whit->
ddl = (*it_gothit)->getDriftDistLeft();
1284 whit->
ddr = (*it_gothit)->getDriftDistRight();
1285 whit->
erddl = (*it_gothit)->getErrDriftDistLeft();
1286 whit->
erddr = (*it_gothit)->getErrDriftDistRight();
1287 whit->
pChiSq = (*it_gothit)->getChisqAdd();
1288 whit->
lr = (*it_gothit)->getFlagLR();
1289 whit->
stat = (*it_gothit)->getStat();
1290 mdcid = (*it_gothit)->getMdcId();
1294 int wid = w0id + localwid;
1296 <<
"lr from PR: "<<whit->
lr
1297 <<
" layerId = " << layid
1298 <<
" wireId = " << localwid
1308 whit->
tdc = (*it_gothit)->getTdc();
1309 whit->
adc= (*it_gothit)->getAdc();
1310 rectrk->
hitcol.push_back(whit);
1311 mhit_mgr->push_back(*whit);
1313 mtrk_mgr->push_back(*rectrk);
1314 mtrkadd_mgr->push_back(*trkadd);
1322 m_trkdelx = trkx1 - trkx2;
1323 m_trkdely = trky1 - trky2;
1324 m_trkdelz = trkz1 - trkz2;
1325 m_trkdelthe = trkthe1 + trkthe2;
1326 m_trkdelphi = trkphi1- trkphi2;
1327 m_trkdelp = trkp1 - trkp2;
1328 StatusCode sc4 = m_nt4->write();
1329 if( sc4.isFailure() ) cout<<
"Ntuple4 filling failed!"<<endl;
1332 if(
debug_ == 4) { std::cout<<
"before refit,ntrk,nhits,nadd..."<<mtrk_mgr->size()
1333 <<
"********"<<mhit_mgr->size()<<
"****"<<mtrkadd_mgr->size()<<endl;
1344 double mdang = 180.0 - csmp3[0].angle(csmp3[1].
unit())*180.0/
M_PI;
1345 double mdphi = 180.0 - fabs(csmphi[0]-csmphi[1])*180.0/
M_PI;
1350 log << MSG::DEBUG <<
"after kalman_fitting(),but in execute...."<<endreq;
1357 SmartDataPtr<RecMdcKalTrackCol> recmdckaltrkCol(eventSvc(),
"/Event/Recon/RecMdcKalTrackCol");
1362 RecMdcKalTrackCol::iterator KalTrk= recmdckaltrkCol->begin();
1364 for(;KalTrk !=recmdckaltrkCol->end();KalTrk++, i_trk++){
1366 bool trkFailed =
false;
1367 for(
int hypo=0; hypo<5; hypo++)
1369 if((*KalTrk)->getStat(0,hypo)==1) {
1370 nFailedTrks[hypo]++;
1374 if(trkFailed) cout<<
" Evt "<<
myEventNo<<
", track "<<i_trk<<
": Kalman filter failed"<<endl;
1377 HelixSegRefVec::iterator iter_hit = gothelixsegs.begin();
1379 int nhitofthistrk=0;
1380 for( ; iter_hit != gothelixsegs.end(); iter_hit++){
1387 iter_hit=gothelixsegs.begin();
1435 SmartDataPtr<Event::CgemMcHitCol> cgemMcHitCol(eventSvc(),
"/Event/MC/CgemMcHitCol");
1437 log << MSG::WARNING <<
"Could not retrieve Cgem MC truth" << endreq;
1438 return StatusCode::FAILURE;
1440 m_evt4=eventHeader->eventNumber();
1441 Event::CgemMcHitCol::iterator iter_truth = cgemMcHitCol->begin();
1443 double dprex,dprey,dprez,dpostx,dposty,dpostz;
1445 for(; iter_truth != cgemMcHitCol->end(); ++iter_truth){
1446 layer = (*iter_truth)->GetLayerID();
1447 dprex = (*iter_truth)->GetPositionXOfPrePoint();
1448 dprey = (*iter_truth)->GetPositionYOfPrePoint();
1449 dprez = (*iter_truth)->GetPositionZOfPrePoint();
1450 dpostx = (*iter_truth)->GetPositionXOfPostPoint();
1451 dposty = (*iter_truth)->GetPositionYOfPostPoint();
1452 dpostz = (*iter_truth)->GetPositionZOfPostPoint();
1460 if((dprex >=0&&dprey>=0)||(dprex <0&&dprey>=0)){
1462 m_dtphi[jj] = acos(dprex/diR);
1465 m_dtpostphi[jj] = acos(dpostx/doR);
1466 m_dtpostz[jj]= dpostz;
1468 if((dprex <0&&dprey<0)||(dprex >=0&&dprey<0)){
1470 m_dtphi[jj] = -acos(dprex/diR);
1473 m_dtpostphi[jj] = -acos(dpostx/doR);
1474 m_dtpostz[jj]= dpostz;
1476 midx=(dprex+dpostx)/2;
1477 midy=(dprey+dposty)/2;
1478 if((midx>=0&&midy>=0)||(midx<0&&midy>=0)){
1480 m_dtphi[jj]=acos(midx/dmR);
1482 if((midx<0&&midy<0)||(midx>=0&&midy<0)){
1484 m_dtphi[jj]=-acos(midx/dmR);
1486 m_dtv[jj] = (m_dtv[jj]+m_dtpostz[jj])/20;
1487 m_tlayer[jj] = layer;
1497 return StatusCode::SUCCESS;
3195 cout<<
"**********************"<<endl;
3196 cout<<
"filter pid type "<<l_mass<<endl;
3201 StatusCode sc = Gaudi::svcLocator()->service(
"MdcGeomSvc", igeomsvc);
3202 if(sc==StatusCode::FAILURE) cout <<
"GeoSvc failing!!!!!!!SC=" << sc << endl;
3205 std::cout<<
"ERROR OCCUR when dynamic_cast in KalFitTrack.cxx !!"<<std::endl;
3209 HepVector pos_old(3,0);
3210 double r0kal_prec(0);
3212 int nhit = track.
HitsMdc().size();
3213 if(
debug_ == 4) cout<<
"filter_fwd..........111 nhit="<<nhit<<endl;
3214 for(
int i=0 ; i < nhit; i++ ) {
3217 if (HitMdc.
chi2()<0)
continue;
3223 int wireid = Wire.
geoID();
3227 Hep3Vector wire = (CLHEP::Hep3Vector)fwd -(CLHEP::Hep3Vector)bck;
3230 work.
pivot((fwd + bck) * .5);
3237 HepPoint3D x0kal = (work.
x(0).z() - bck.z())/ wire.z() * wire + bck;
3240 std::cout<<
" a="<<track.
a()<<std::endl;
3241 std::cout<<
" pivot="<<track.
pivot()<<std::endl;
3242 std::cout<<
" x0kal before sag: "<<x0kal<<std::endl;
3287 double tension = geowire->
Tension();
3289 double zinit(x0kal.z()), lzx(Wire.
lzx());
3291 double A = 47.35E-6/tension;
3292 double Zp = (zinit - bck.z())*lzx/wire.z();
3295 std::cout<<
" sag in filter_fwd_anal: "<<std::endl;
3296 std::cout<<
" x0kal.x(): "<<std::setprecision(10)<<x0kal.x()<<std::endl;
3297 std::cout<<
"zinit: "<<zinit<<
" bck.z(): "<<bck.z()<<std::endl;
3298 std::cout<<
" wire.x()*(zinit-bck.z())/wire.z(): "<<std::setprecision(10)
3299 <<(wire.x()*(zinit-bck.z())/wire.z())<<std::endl;
3300 std::cout<<
"bck.x(): "<<std::setprecision(10)<<bck.x()<<std::endl;
3301 std::cout<<
" wire.x()*(zinit-bck.z())/wire.z() + bck.x(): "<<std::setprecision(10)
3302 <<(wire.x()*(zinit-bck.z())/wire.z() + bck.x())<<std::endl;
3305 result.setX(wire.x()*(zinit-bck.z())/wire.z() + bck.x());
3306 result.setY((A*(Zp-lzx)+wire.y()/lzx)*Zp+bck.y());
3307 result.setZ((A*(2*Zp-lzx)*lzx+wire.y())/wire.z());
3309 wire.setX(wire.x()/wire.z());
3310 wire.setY(result.z());
3312 x0kal.setX(result.x());
3313 x0kal.setY(result.y());
3318 double r0kal = x0kal.perp();
3323 if(i==0) track.
pivot(x0kal);
3327 if (nhits_read == 1) {
3338 std::cout<<
" x0kal after sag: "<<x0kal<<std::endl;
3339 std::cout<<
" a="<<track.
a()<<std::endl;
3340 std::cout<<
" Ea="<<track.
Ea()<<std::endl;
3343 double dtracknew = 0.;
3347 if(fabs(track.
kappa())>0&&fabs(track.
kappa())<1000.0&&fabs(track.
tanl())<7.02) {
3348 Hep3Vector meas = track.
momentum(0).cross(wire).unit();
3349 double diff_chi2 = track.
chiSq();
3350 Hep3Vector IP(0,0,0);
3356 for(
unsigned k=i+1 ; k < nhit; k++ )
3361 double dchi2 = -1.0;
3416 m_residest = dtrack-dtdc;
3417 m_residnew = dtracknew -dtdc;
3420 m_anal_dr = worktemp.
a()[0];
3421 m_anal_phi0 = worktemp.
a()[1];
3422 m_anal_kappa = worktemp.
a()[2];
3423 m_anal_dz = worktemp.
a()[3];
3424 m_anal_tanl = worktemp.
a()[4];
3425 m_anal_ea_dr = worktemp.
Ea()[0][0];
3426 m_anal_ea_phi0 = worktemp.
Ea()[1][1];
3427 m_anal_ea_kappa = worktemp.
Ea()[2][2];
3428 m_anal_ea_dz = worktemp.
Ea()[3][3];
3429 m_anal_ea_tanl = worktemp.
Ea()[4][4];
3430 StatusCode sc5 = m_nt5->write();
3431 if( sc5.isFailure() ) cout<<
"Ntuple5 filling failed!"<<endl;
3439 m_dchi2_hit_e[m_hit_no] = dchi2;
3440 m_residest_hit_e[m_hit_no] = dtrack-dtdc;
3441 m_residnew_hit_e[m_hit_no] = dtracknew -dtdc;
3444 m_anal_dr_hit_e[m_hit_no] = worktemp.
a()[0];
3445 m_anal_phi0_hit_e[m_hit_no] = worktemp.
a()[1];
3446 m_anal_kappa_hit_e[m_hit_no] = worktemp.
a()[2];
3447 m_anal_dz_hit_e[m_hit_no] = worktemp.
a()[3];
3448 m_anal_tanl_hit_e[m_hit_no] = worktemp.
a()[4];
3449 m_anal_ea_dr_hit_e[m_hit_no] = worktemp.
Ea()[0][0];
3450 m_anal_ea_phi0_hit_e[m_hit_no] = worktemp.
Ea()[1][1];
3451 m_anal_ea_kappa_hit_e[m_hit_no] = worktemp.
Ea()[2][2];
3452 m_anal_ea_dz_hit_e[m_hit_no] = worktemp.
Ea()[3][3];
3453 m_anal_ea_tanl_hit_e[m_hit_no] = worktemp.
Ea()[4][4];
3457 else if(l_mass == 1){
3458 m_dchi2_hit_mu[m_hit_no] = dchi2;
3459 m_residest_hit_mu[m_hit_no] = dtrack-dtdc;
3460 m_residnew_hit_mu[m_hit_no] = dtracknew -dtdc;
3463 m_anal_dr_hit_mu[m_hit_no] = worktemp.
a()[0];
3464 m_anal_phi0_hit_mu[m_hit_no] = worktemp.
a()[1];
3465 m_anal_kappa_hit_mu[m_hit_no] = worktemp.
a()[2];
3466 m_anal_dz_hit_mu[m_hit_no] = worktemp.
a()[3];
3467 m_anal_tanl_hit_mu[m_hit_no] = worktemp.
a()[4];
3468 m_anal_ea_dr_hit_mu[m_hit_no] = worktemp.
Ea()[0][0];
3469 m_anal_ea_phi0_hit_mu[m_hit_no] = worktemp.
Ea()[1][1];
3470 m_anal_ea_kappa_hit_mu[m_hit_no] = worktemp.
Ea()[2][2];
3471 m_anal_ea_dz_hit_mu[m_hit_no] = worktemp.
Ea()[3][3];
3472 m_anal_ea_tanl_hit_mu[m_hit_no] = worktemp.
Ea()[4][4];
3477 else if(l_mass == 2){
3478 m_dchi2_hit_pi[m_hit_no] = dchi2;
3479 m_residest_hit_pi[m_hit_no] = dtrack-dtdc;
3480 m_residnew_hit_pi[m_hit_no] = dtracknew -dtdc;
3483 m_anal_dr_hit_pi[m_hit_no] = worktemp.
a()[0];
3484 m_anal_phi0_hit_pi[m_hit_no] = worktemp.
a()[1];
3485 m_anal_kappa_hit_pi[m_hit_no] = worktemp.
a()[2];
3486 m_anal_dz_hit_pi[m_hit_no] = worktemp.
a()[3];
3487 m_anal_tanl_hit_pi[m_hit_no] = worktemp.
a()[4];
3488 m_anal_ea_dr_hit_pi[m_hit_no] = worktemp.
Ea()[0][0];
3489 m_anal_ea_phi0_hit_pi[m_hit_no] = worktemp.
Ea()[1][1];
3490 m_anal_ea_kappa_hit_pi[m_hit_no] = worktemp.
Ea()[2][2];
3491 m_anal_ea_dz_hit_pi[m_hit_no] = worktemp.
Ea()[3][3];
3492 m_anal_ea_tanl_hit_pi[m_hit_no] = worktemp.
Ea()[4][4];
3497 else if(l_mass == 3){
3498 m_dchi2_hit_k[m_hit_no] = dchi2;
3499 m_residest_hit_k[m_hit_no] = dtrack-dtdc;
3500 m_residnew_hit_k[m_hit_no] = dtracknew -dtdc;
3503 m_anal_dr_hit_k[m_hit_no] = worktemp.
a()[0];
3504 m_anal_phi0_hit_k[m_hit_no] = worktemp.
a()[1];
3505 m_anal_kappa_hit_k[m_hit_no] = worktemp.
a()[2];
3506 m_anal_dz_hit_k[m_hit_no] = worktemp.
a()[3];
3507 m_anal_tanl_hit_k[m_hit_no] = worktemp.
a()[4];
3508 m_anal_ea_dr_hit_k[m_hit_no] = worktemp.
Ea()[0][0];
3509 m_anal_ea_phi0_hit_k[m_hit_no] = worktemp.
Ea()[1][1];
3510 m_anal_ea_kappa_hit_k[m_hit_no] = worktemp.
Ea()[2][2];
3511 m_anal_ea_dz_hit_k[m_hit_no] = worktemp.
Ea()[3][3];
3512 m_anal_ea_tanl_hit_k[m_hit_no] = worktemp.
Ea()[4][4];
3517 else if(l_mass == 4){
3518 m_dchi2_hit_p[m_hit_no] = dchi2;
3519 m_residest_hit_p[m_hit_no] = dtrack-dtdc;
3520 m_residnew_hit_p[m_hit_no] = dtracknew -dtdc;
3523 m_anal_dr_hit_p[m_hit_no] = worktemp.
a()[0];
3524 m_anal_phi0_hit_p[m_hit_no] = worktemp.
a()[1];
3525 m_anal_kappa_hit_p[m_hit_no] = worktemp.
a()[2];
3526 m_anal_dz_hit_p[m_hit_no] = worktemp.
a()[3];
3527 m_anal_tanl_hit_p[m_hit_no] = worktemp.
a()[4];
3528 m_anal_ea_dr_hit_p[m_hit_no] = worktemp.
Ea()[0][0];
3529 m_anal_ea_phi0_hit_p[m_hit_no] = worktemp.
Ea()[1][1];
3530 m_anal_ea_kappa_hit_p[m_hit_no] = worktemp.
Ea()[2][2];
3531 m_anal_ea_dz_hit_p[m_hit_no] = worktemp.
Ea()[3][3];
3532 m_anal_ea_tanl_hit_p[m_hit_no] = worktemp.
Ea()[4][4];
3543 diff_chi2 = chi2 - diff_chi2;
3544 HitMdc.
chi2(diff_chi2);
3854 MsgStream log(
msgSvc(), name());
3855 double Pt_threshold(0.3);
3856 Hep3Vector IP(0,0,0);
3862 if ( !&whMgr )
return;
3865 int ntrk = mdcMgr->size();
3866 double* rPt =
new double[ntrk];
3867 int* rOM =
new int[ntrk];
3868 unsigned int* order =
new unsigned int[ntrk];
3869 unsigned int* rCont =
new unsigned int[ntrk];
3870 unsigned int* rGen =
new unsigned int[ntrk];
3872 if(
testOutput) cout<<
"nMdcTrk = "<<ntrk<<endl;
3875 for(vector<MdcRec_trk>::iterator it = mdcMgr->begin(),
3876 end = mdcMgr->end(); it != end; it++) {
3880 rPt[index] = 1 / fabs(it->helix[2]);
3881 if(
debug_ == 4) cout<<
"rPt...[ "<<index<<
" ]...."<< rPt[index] <<endl;
3882 if(rPt[index] < 0) rPt[index] =
DBL_MAX;
3884 std::vector<MdcRec_wirhit*> pt = it->hitcol ;
3885 if(
debug_ == 4) cout<<
"ppt size: "<< pt.size()<<endl;
3887 for (vector<MdcRec_wirhit*>::iterator ii = pt.end()-1;
3888 ii !=pt.begin()-1; ii--) {
3889 int lyr((*ii)->geo->Lyr()->Id());
3890 if (outermost < lyr) outermost = lyr;
3891 if(
debug_ == 4) cout<<
"outmost: "<<outermost<<
" lyr: "<<lyr<<endl;
3893 rOM[index] = outermost;
3894 order[index] = index;
3899 for (
int j, k = ntrk - 1; k >= 0; k = j){
3901 for(
int i = 1; i <= k; i++)
3902 if(rPt[i - 1] < rPt[i]){
3904 std::swap(order[i], order[j]);
3905 std::swap(rPt[i], rPt[j]);
3906 std::swap(rOM[i], rOM[j]);
3907 std::swap(rCont[i], rCont[j]);
3908 std::swap(rGen[i], rGen[j]);
3916 DataObject *aReconEvent;
3917 eventSvc()->findObject(
"/Event/Recon",aReconEvent);
3921 StatusCode sc = eventSvc()->registerObject(
"/Event/Recon",recevt );
3922 if(sc!=StatusCode::SUCCESS) {
3923 log << MSG::FATAL <<
"Could not register ReconEvent" <<endreq;
3931 log << MSG::INFO <<
"beginning to make RecMdcKalTrackCol" <<endreq;
3933 for(
int l = 0; l < ntrk; l++) {
3942 MdcRec_trk& TrasanTRK = *(mdcMgr->begin() + order[l]);
3946 int trasqual = TrasanTRK_add.
quality;
3947 if(
debug_ == 4) cout<<
"kalman_fitting>trasqual... "<<trasqual<<endl;
3948 if (trasqual)
continue;
3952 cout <<
"******* KalFit NUMBER : " << newcount << std::endl;
3956 if ((TrasanTRK_add.
decision & 32) == 32 ||
3957 (TrasanTRK_add.
decision & 64) == 64) type = 1;
3962 TrasanTRK.
pivot[2]);
3965 for(
int i = 0; i < 5; i++)
3966 a[i] = TrasanTRK.
helix[i];
3969 for(
int i = 0, k = 0; i < 5; i++) {
3970 for(
int j = 0; j <= i; j++) {
3972 ea[j][i] = ea[i][j];
3975 double fiTerm = TrasanTRK.
fiTerm;
3981 double pT_trk = fabs(track_lead.
pt());
3984 int inlyr(999), outlyr(-1);
3985 int* rStat =
new int[43];
3986 for(
int irStat=0;irStat<43;++irStat)rStat[irStat]=0;
3987 std::vector<MdcRec_wirhit*> pt=TrasanTRK.
hitcol;
3989 if(
debug_ == 4) cout<<
"*********Pt size****"<< pt.size()<<endl;
3992 for (vector<MdcRec_wirhit*>::iterator ii = pt.end()-1;
3993 ii != pt.begin()-1; ii--) {
3994 Num[(*ii)->geo->Lyr()->Id()]++;
3997 for (vector<MdcRec_wirhit*>::iterator ii = pt.end()-1;
3998 ii != pt.begin()-1; ii--) {
4001 if (pT_trk>0.12 && Num[(*ii)->geo->Lyr()->Id()]>3) {
4003 cout <<
"WARNING: I found " << Num[(*ii)->geo->Lyr()->Id()]
4004 <<
" hits in the layer "
4005 << (*ii)->geo->Lyr()->Id() << std::endl;
4025 double dist[2] = {rechit.
ddl, rechit.
ddr};
4026 double erdist[2] = {rechit.
erddl, rechit.
erddr};
4031 if (rechit.
lr==2 || rechit.
lr==0) lr_decision=-1;
4033 else if (rechit.
lr==1) lr_decision=1;
4036 int ind(geo->
Lyr()->
Id());
4040 lr_decision, rechit.
tdc,
4045 if (inlyr>ind) inlyr = ind;
4046 if (outlyr<ind) outlyr = ind;
4049 cout <<
"**** NUMBER OF Mdc HITS (TRASAN) = " << hit_asso << std::endl;
4052 int empty_between(0), empty(0);
4053 for (
int i= inlyr; i <= outlyr; i++)
4054 if (!rStat[i]) empty_between++;
4055 empty = empty_between+inlyr+(42-outlyr);
4062 for(std::vector<KalFitHitMdc>::iterator it_hit = track_lead.
HitsMdc().begin(); it_hit!=track_lead.
HitsMdc().end(); it_hit++){
4068 track_lead.
type(type);
4070 unsigned int nhit = track_lead.
HitsMdc().size();
4071 if (!nhit &&
debug_ == 4) {
4072 cout <<
" ATTENTION TRACK WITH ONLY HITS " << nhit << std::endl;
4077 double KalFitst(0), KalFitax(0), KalFitschi2(0);
4080 Hep3Vector outer_pivot(track_lead.
x(fiTerm));
4081 track_lead.
pivot(outer_pivot);
4087 HepSymMatrix Eakal(5,0);
4090 double costheta = track_lead.
a()[4] / sqrt(1.0 + track_lead.
a()[4]*track_lead.
a()[4]);
4100 cout <<
"from Mdc Pattern Recognition: " << std::endl;
4106 cout <<
" dr = " << work.
a()[0]
4107 <<
", Er_dr = " << sqrt(work.
Ea()[0][0]) << std::endl;
4108 cout <<
" phi0 = " << work.
a()[1]
4109 <<
", Er_phi0 = " << sqrt(work.
Ea()[1][1]) << std::endl;
4110 cout <<
" PT = " << 1/work.
a()[2]
4111 <<
", Er_kappa = " << sqrt(work.
Ea()[2][2]) << std::endl;
4112 cout <<
" dz = " << work.
a()[3]
4113 <<
", Er_dz = " << sqrt(work.
Ea()[3][3]) << std::endl;
4114 cout <<
" tanl = " << work.
a()[4]
4115 <<
", Er_tanl = " << sqrt(work.
Ea()[4][4]) << std::endl;
4119 cout<<
"----------------------"<<endl;
4120 cout<<
"track "<<l<<
", pid = "<<
lead_<<
": "<<endl;
4137 cout <<
" Mdc FIRST KALMAN FIT " << std::endl;
4142 cout <<
" dr = " << work.
a()[0]
4143 <<
", Er_dr = " << sqrt(work.
Ea()[0][0]) << std::endl;
4144 cout <<
" phi0 = " << work.
a()[1]
4145 <<
", Er_phi0 = " << sqrt(work.
Ea()[1][1]) << std::endl;
4146 cout <<
" PT = " << 1/work.
a()[2]
4147 <<
", Er_kappa = " << sqrt(work.
Ea()[2][2]) << std::endl;
4148 cout <<
" dz = " << work.
a()[3]
4149 <<
", Er_dz = " << sqrt(work.
Ea()[3][3]) << std::endl;
4150 cout <<
" tanl = " << work.
a()[4]
4151 <<
", Er_tanl = " << sqrt(work.
Ea()[4][4]) << std::endl;
4158 complete_track(TrasanTRK, TrasanTRK_add, track_lead, kaltrk,kalcol,segcol,1);
4162 IDataProviderSvc* eventSvc = NULL;
4163 Gaudi::svcLocator()->service(
"EventDataSvc", eventSvc);
4165 log << MSG::INFO <<
"makeTds:event Svc has been found" << endreq;
4167 log << MSG::FATAL <<
"makeTds:Could not find eventSvc" << endreq;
4171 IDataManagerSvc *dataManSvc;
4172 dataManSvc=
dynamic_cast<IDataManagerSvc*
>(eventSvc);
4173 DataObject *aKalTrackCol;
4174 eventSvc->findObject(
"/Event/Recon/RecMdcKalTrackCol",aKalTrackCol);
4175 if(aKalTrackCol != NULL) {
4176 dataManSvc->clearSubTree(
"/Event/Recon/RecMdcKalTrackCol");
4177 eventSvc->unregisterObject(
"/Event/Recon/RecMdcKalTrackCol");
4179 kalsc = eventSvc->registerObject(
"/Event/Recon/RecMdcKalTrackCol", kalcol);
4180 if( kalsc.isFailure() ) {
4181 log << MSG::FATAL <<
"Could not register RecMdcKalTrack" << endreq;
4184 log << MSG::INFO <<
"RecMdcKalTrackCol registered successfully!" <<endreq;
4187 DataObject *aRecKalSegEvent;
4188 eventSvc->findObject(
"/Event/Recon/RecMdcKalHelixSegCol", aRecKalSegEvent);
4189 if(aRecKalSegEvent!=NULL) {
4191 segsc = eventSvc->unregisterObject(
"/Event/Recon/RecMdcKalHelixSegCol");
4192 if(segsc != StatusCode::SUCCESS) {
4193 log << MSG::FATAL <<
"Could not unregister RecMdcKalHelixSegCol collection" << endreq;
4197 segsc = eventSvc->registerObject(
"/Event/Recon/RecMdcKalHelixSegCol", segcol);
4198 if( segsc.isFailure() ) {
4199 log << MSG::FATAL <<
"Could not register RecMdcKalHelixSeg" << endreq;
4202 log << MSG::INFO <<
"RecMdcKalHelixSegCol registered successfully!" <<endreq;
4204 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.);
4205 double r1(0.),r2(0.),kap1(999.),kap2(999.),tanl1(0),tanl2(0.);
4206 double px1(0.),px2(0.),py1(0.),py2(0.),pz1(0.),pz2(0.),charge1(0.),charge2(0.);
4209 SmartDataPtr<RecMdcKalTrackCol> kaltrkCol(eventSvc,
"/Event/Recon/RecMdcKalTrackCol");
4211 log << MSG::FATAL <<
"Could not find RecMdcKalTrackCol" << endreq;
4214 log << MSG::INFO <<
"Begin to check RecMdcKalTrackCol"<<endreq;
4215 RecMdcKalTrackCol::iterator iter_trk = kaltrkCol->begin();
4216 for(
int jj=1; iter_trk != kaltrkCol->end(); iter_trk++,jj++) {
4217 log << MSG::DEBUG <<
"retrieved MDC Kalmantrack:"
4218 <<
"Track Id: " << (*iter_trk)->getTrackId()
4219 <<
" Mass of the fit: "<< (*iter_trk)->getMass(2)<< endreq
4220 <<
" Length of the track: "<< (*iter_trk)->getLength(2)
4221 <<
" Tof of the track: "<< (*iter_trk)->getTof(2) << endreq
4222 <<
" Chisq of the fit: "<< (*iter_trk)->getChisq(0,2)
4223 <<
" "<< (*iter_trk)->getChisq(1,2) << endreq
4224 <<
"Ndf of the fit: "<< (*iter_trk)->getNdf(0,2)
4225 <<
" "<< (*iter_trk)->getNdf(1,2) << endreq
4226 <<
"Kappa " << (*iter_trk)->getZHelix()[2]
4228 for(
int i = 0; i<43; i++) {
4229 log << MSG::DEBUG <<
"retrieved pathl["<<i<<
"]= "
4230 << (*iter_trk)->getPathl(i) <<endreq;
4233 m_trackid = (*iter_trk)->getTrackId();
4235 for(
int jj =0, iii=0; jj<5; jj++) {
4237 m_length[jj] = (*iter_trk)->getLength(jj);
4238 m_tof[jj] = (*iter_trk)->getTof(jj);
4239 m_nhits[jj] = (*iter_trk)->getNhits(jj);
4240 m_zhelix[jj] = (*iter_trk)->getZHelix()[jj];
4241 m_zhelixe[jj] = (*iter_trk)->getZHelixE()[jj];
4242 m_zhelixmu[jj] = (*iter_trk)->getZHelixMu()[jj];
4243 m_zhelixk[jj] = (*iter_trk)->getZHelixK()[jj];
4244 m_zhelixp[jj] = (*iter_trk)->getZHelixP()[jj];
4245 m_fhelix[jj] = (*iter_trk)->getFHelix()[jj];
4246 m_fhelixe[jj] = (*iter_trk)->getFHelixE()[jj];
4247 m_fhelixmu[jj] = (*iter_trk)->getFHelixMu()[jj];
4248 m_fhelixk[jj] = (*iter_trk)->getFHelixK()[jj];
4249 m_fhelixp[jj] = (*iter_trk)->getFHelixP()[jj];
4250 m_lhelix[jj] = (*iter_trk)->getLHelix()[jj];
4251 m_lhelixe[jj] = (*iter_trk)->getLHelixE()[jj];
4252 m_lhelixmu[jj] = (*iter_trk)->getLHelixMu()[jj];
4253 m_lhelixk[jj] = (*iter_trk)->getLHelixK()[jj];
4254 m_lhelixp[jj] = (*iter_trk)->getLHelixP()[jj];
4257 for(
int kk=0; kk<=jj; kk++,iii++) {
4258 m_zerror[iii] = (*iter_trk)->getZError()[jj][kk];
4259 m_zerrore[iii] = (*iter_trk)->getZErrorE()[jj][kk];
4260 m_zerrormu[iii] = (*iter_trk)->getZErrorMu()[jj][kk];
4261 m_zerrork[iii] = (*iter_trk)->getZErrorK()[jj][kk];
4262 m_zerrorp[iii] = (*iter_trk)->getZErrorP()[jj][kk];
4263 m_ferror[iii] = (*iter_trk)->getFError()[jj][kk];
4264 m_ferrore[iii] = (*iter_trk)->getFErrorE()[jj][kk];
4265 m_ferrormu[iii] = (*iter_trk)->getFErrorMu()[jj][kk];
4266 m_ferrork[iii] = (*iter_trk)->getFErrorK()[jj][kk];
4267 m_ferrorp[iii] = (*iter_trk)->getFErrorP()[jj][kk];
4268 m_lerror[iii] = (*iter_trk)->getLError()[jj][kk];
4269 m_lerrore[iii] = (*iter_trk)->getLErrorE()[jj][kk];
4270 m_lerrormu[iii] = (*iter_trk)->getLErrorMu()[jj][kk];
4271 m_lerrork[iii] = (*iter_trk)->getLErrorK()[jj][kk];
4272 m_lerrorp[iii] = (*iter_trk)->getLErrorP()[jj][kk];
4312 m_chisq[0][0] = (*iter_trk)->getChisq(0,0);
4313 m_chisq[1][0] = (*iter_trk)->getChisq(0,1);
4314 m_chisq[2][0] = (*iter_trk)->getChisq(0,2);
4315 m_chisq[3][0] = (*iter_trk)->getChisq(0,3);
4316 m_chisq[4][0] = (*iter_trk)->getChisq(0,4);
4317 m_chisq[0][1] = (*iter_trk)->getChisq(1,0);
4318 m_chisq[1][1] = (*iter_trk)->getChisq(1,1);
4319 m_chisq[2][1] = (*iter_trk)->getChisq(1,2);
4320 m_chisq[3][1] = (*iter_trk)->getChisq(1,3);
4321 m_chisq[4][1] = (*iter_trk)->getChisq(1,4);
4323 m_ndf[0][0] = (*iter_trk)->getNdf(0,0);
4324 m_ndf[1][0] = (*iter_trk)->getNdf(0,1);
4325 m_ndf[2][0] = (*iter_trk)->getNdf(0,2);
4326 m_ndf[3][0] = (*iter_trk)->getNdf(0,3);
4327 m_ndf[4][0] = (*iter_trk)->getNdf(0,4);
4328 m_ndf[0][1] = (*iter_trk)->getNdf(1,0);
4329 m_ndf[1][1] = (*iter_trk)->getNdf(1,1);
4330 m_ndf[2][1] = (*iter_trk)->getNdf(1,2);
4331 m_ndf[3][1] = (*iter_trk)->getNdf(1,3);
4332 m_ndf[4][1] = (*iter_trk)->getNdf(1,4);
4334 m_stat[0][0] = (*iter_trk)->getStat(0,0);
4335 m_stat[1][0] = (*iter_trk)->getStat(0,1);
4336 m_stat[2][0] = (*iter_trk)->getStat(0,2);
4337 m_stat[3][0] = (*iter_trk)->getStat(0,3);
4338 m_stat[4][0] = (*iter_trk)->getStat(0,4);
4339 m_stat[0][1] = (*iter_trk)->getStat(1,0);
4340 m_stat[1][1] = (*iter_trk)->getStat(1,1);
4341 m_stat[2][1] = (*iter_trk)->getStat(1,2);
4342 m_stat[3][1] = (*iter_trk)->getStat(1,3);
4343 m_stat[4][1] = (*iter_trk)->getStat(1,4);
4345 m_fptot = sqrt(1+pow(m_fhelix[4],2))/m_fhelix[2];
4346 m_fptote = sqrt(1+pow(m_fhelixe[4],2))/m_fhelixe[2];
4347 m_fptotmu = sqrt(1+pow(m_fhelixmu[4],2))/m_fhelixmu[2];
4348 m_fptotk = sqrt(1+pow(m_fhelixk[4],2))/m_fhelixk[2];
4349 m_fptotp = sqrt(1+pow(m_fhelixp[4],2))/m_fhelixp[2];
4351 m_lptot = sqrt(1+pow(m_lhelix[4],2))/m_lhelix[2];
4352 m_lptote = sqrt(1+pow(m_lhelixe[4],2))/m_lhelixe[2];
4353 m_lptotmu = sqrt(1+pow(m_lhelixmu[4],2))/m_lhelixmu[2];
4354 m_lptotk = sqrt(1+pow(m_lhelixk[4],2))/m_lhelixk[2];
4355 m_lptotp = sqrt(1+pow(m_lhelixp[4],2))/m_lhelixp[2];
4357 m_lpt = 1/m_lhelix[2];
4358 m_lpte = 1/m_lhelixe[2];
4359 m_lptmu = 1/m_lhelixmu[2];
4360 m_lptk = 1/m_lhelixk[2];
4361 m_lptp = 1/m_lhelixp[2];
4363 m_fpt = 1/m_fhelix[2];
4364 m_fpte = 1/m_fhelixe[2];
4365 m_fptmu = 1/m_fhelixmu[2];
4366 m_fptk = 1/m_fhelixk[2];
4367 m_fptp = 1/m_fhelixp[2];
4370 std::cout<<
" "<<std::endl;
4371 std::cout<<
"in file Kalman_fitting_anal ,the m_fpt is .."<<m_fpt<<std::endl;
4372 std::cout<<
"in file Kalman_fitting_anal ,the m_fpte is .."<<m_fpte<<std::endl;
4373 std::cout<<
"in file Kalman_fitting_anal ,the m_fptmu is .."<<m_fptmu<<std::endl;
4374 std::cout<<
"in file Kalman_fitting_anal ,the m_fptk is .."<<m_fptk<<std::endl;
4375 std::cout<<
"in file Kalman_fitting_anal ,the m_fptp is .."<<m_fptp<<std::endl;
4378 m_zpt = 1/m_zhelix[2];
4379 m_zpte = 1/m_zhelixe[2];
4380 m_zptmu = 1/m_zhelixmu[2];
4381 m_zptk = 1/m_zhelixk[2];
4382 m_zptp = 1/m_zhelixp[2];
4385 std::cout<<
"in file Kalman_fitting_anal ,the m_zpt is .."<<m_zpt<<std::endl;
4386 std::cout<<
"in file Kalman_fitting_anal ,the m_zpte is .."<<m_zpte<<std::endl;
4387 std::cout<<
"in file Kalman_fitting_anal ,the m_zptmu is .."<<m_zptmu<<std::endl;
4388 std::cout<<
"in file Kalman_fitting_anal ,the m_zptk is .."<<m_zptk<<std::endl;
4389 std::cout<<
"in file Kalman_fitting_anal ,the m_zptp is .."<<m_zptp<<std::endl;
4391 m_zptot = sqrt(1+pow(m_zhelix[4],2))/m_zhelix[2];
4392 m_zptote = sqrt(1+pow(m_zhelixe[4],2))/m_zhelixe[2];
4393 m_zptotmu = sqrt(1+pow(m_zhelixmu[4],2))/m_zhelixmu[2];
4394 m_zptotk = sqrt(1+pow(m_zhelixk[4],2))/m_zhelixk[2];
4395 m_zptotp = sqrt(1+pow(m_zhelixp[4],2))/m_zhelixp[2];
4398 std::cout<<
"in file Kalman_fitting_anal ,the m_zptot is .."<<m_zptot<<std::endl;
4399 std::cout<<
"in file Kalman_fitting_anal ,the m_zptote is .."<<m_zptote<<std::endl;
4400 std::cout<<
"in file Kalman_fitting_anal ,the m_zptotmu is .."<<m_zptotmu<<std::endl;
4401 std::cout<<
"in file Kalman_fitting_anal ,the m_zptotk is .."<<m_zptotk<<std::endl;
4402 std::cout<<
"in file Kalman_fitting_anal ,the m_zptotp is .."<<m_zptotp<<std::endl;
4406 m_zsigp = sqrt(pow((m_zptot/m_zhelix[2]),2)*m_zerror[5]+
4407 pow((m_zhelix[4]/m_zptot),2)*pow((1/m_zhelix[2]),4)*m_zerror[14]-
4408 2*m_zhelix[4]*m_zerror[12]*pow((1/m_zhelix[2]),3));
4409 m_zsigpe = sqrt(pow((m_zptote/m_zhelixe[2]),2)*m_zerrore[5]+
4410 pow((m_zhelixe[4]/m_zptote),2)*pow((1/m_zhelixe[2]),4)*m_zerrore[14]-
4411 2*m_zhelixe[4]*m_zerrore[12]*pow((1/m_zhelixe[2]),3));
4412 m_zsigpmu = sqrt(pow((m_zptotmu/m_zhelixmu[2]),2)*m_zerrormu[5]+
4413 pow((m_zhelixmu[4]/m_zptotmu),2)*pow((1/m_zhelixmu[2]),4)*m_zerrormu[14]-
4414 2*m_zhelixmu[4]*m_zerrormu[12]*pow((1/m_zhelixmu[2]),3));
4415 m_zsigpk = sqrt(pow((m_zptotk/m_zhelixk[2]),2)*m_zerrork[5]+
4416 pow((m_zhelixk[4]/m_zptotk),2)*pow((1/m_zhelixk[2]),4)*m_zerrork[14]-
4417 2*m_zhelixk[4]*m_zerrork[12]*pow((1/m_zhelixk[2]),3));
4418 m_zsigpp = sqrt(pow((m_zptotp/m_zhelixp[2]),2)*m_zerrorp[5]+
4419 pow((m_zhelixp[4]/m_zptotp),2)*pow((1/m_zhelixp[2]),4)*m_zerrorp[14]-
4420 2*m_zhelixp[4]*m_zerrorp[12]*pow((1/m_zhelixp[2]),3));
4423 StatusCode sc1 = m_nt1->write();
4424 if( sc1.isFailure() ) cout<<
"Ntuple1 filling failed!"<<endl;
4429 phi1 = (*iter_trk)->getFFi0();
4430 r1 = (*iter_trk)->getFDr();
4431 z1 = (*iter_trk)->getFDz();
4432 kap1 = (*iter_trk)->getFCpa();
4433 tanl1 = (*iter_trk)->getFTanl();
4434 charge1 = kap1/fabs(kap1);
4437 p1 = sqrt(1+tanl1*tanl1)/kap1;
4438 the1 =
M_PI/2-atan(tanl1);
4441 pz1= tanl1/fabs(kap1);
4443 }
else if(jj == 2) {
4444 phi2 = (*iter_trk)->getFFi0();
4445 r2 = (*iter_trk)->getFDr();
4446 z2 = (*iter_trk)->getFDz();
4447 kap2 = (*iter_trk)->getFCpa();
4448 tanl2 = (*iter_trk)->getFTanl();
4449 charge2 = kap2/fabs(kap2);
4452 p2 = sqrt(1+tanl2*tanl2)/kap1;
4453 the2 =
M_PI/2-atan(tanl2);
4456 pz2= tanl2/fabs(kap2);
4464 m_delthe = the1 + the2;
4467 m_delpx = charge1*fabs(px1) + charge2*fabs(px2);
4468 m_delpy = charge1*fabs(py1) + charge2*fabs(py2);
4469 m_delpz = charge1*fabs(pz1) + charge2*fabs(pz2);
4471 StatusCode sc2 = m_nt2->write();
4472 if( sc2.isFailure() ) cout<<
"Ntuple2 filling failed!"<<endl;
4481 cout <<
"Kalfitting finished " << std::endl;
4488 MsgStream log(
msgSvc(), name());
4489 double Pt_threshold(0.3);
4490 Hep3Vector IP(0,0,0);
4497 if ( !&whMgr )
return;
4500 int ntrk = mdcMgr->size();
4501 double* rPt =
new double[ntrk];
4502 int* rOM =
new int[ntrk];
4503 unsigned int* order =
new unsigned int[ntrk];
4504 unsigned int* rCont =
new unsigned int[ntrk];
4505 unsigned int* rGen =
new unsigned int[ntrk];
4508 for(vector<MdcRec_trk>::iterator it = mdcMgr->begin(),
4509 end = mdcMgr->end(); it != end; it++) {
4513 rPt[index] = 1 / fabs(it->helix[2]);
4514 if(
debug_ == 4) cout<<
"rPt...[ "<<index<<
" ]...."<< rPt[index] <<endl;
4515 if(rPt[index] < 0) rPt[index] =
DBL_MAX;
4517 std::vector<MdcRec_wirhit*> pt = it->hitcol ;
4518 if(
debug_ == 4) cout<<
"ppt size: "<< pt.size()<<endl;
4520 for (vector<MdcRec_wirhit*>::iterator ii = pt.end()-1;
4521 ii !=pt.begin()-1; ii--) {
4522 int lyr((*ii)->geo->Lyr()->Id());
4523 if (outermost < lyr) outermost = lyr;
4524 if(
debug_ == 4) cout<<
"outmost: "<<outermost<<
" lyr: "<<lyr<<endl;
4526 rOM[index] = outermost;
4527 order[index] = index;
4532 for (
int j, k = ntrk - 1; k >= 0; k = j){
4534 for(
int i = 1; i <= k; i++)
4535 if(rPt[i - 1] < rPt[i]){
4537 std::swap(order[i], order[j]);
4538 std::swap(rPt[i], rPt[j]);
4539 std::swap(rOM[i], rOM[j]);
4540 std::swap(rCont[i], rCont[j]);
4541 std::swap(rGen[i], rGen[j]);
4548 DataObject *aReconEvent;
4549 eventSvc()->findObject(
"/Event/Recon",aReconEvent);
4553 StatusCode sc = eventSvc()->registerObject(
"/Event/Recon",recevt );
4554 if(sc!=StatusCode::SUCCESS) {
4555 log << MSG::FATAL <<
"Could not register ReconEvent" <<endreq;
4563 log << MSG::INFO <<
"beginning to make RecMdcKalTrackCol" <<endreq;
4566 for(
int l = 0; l < ntrk; l++) {
4568 MdcRec_trk& TrasanTRK = *(mdcMgr->begin() + order[l]);
4572 int trasqual = TrasanTRK_add.
quality;
4573 if(
debug_ == 4) cout<<
"kalman_fitting>trasqual... "<<trasqual<<endl;
4574 if (trasqual)
continue;
4578 cout <<
"******* KalFit NUMBER : " << newcount << std::endl;
4582 if ((TrasanTRK_add.
decision & 32) == 32 ||
4583 (TrasanTRK_add.
decision & 64) == 64) type = 1;
4588 TrasanTRK.
pivot[2]);
4591 for(
int i = 0; i < 5; i++)
4592 a[i] = TrasanTRK.
helix[i];
4595 for(
int i = 0, k = 0; i < 5; i++) {
4596 for(
int j = 0; j <= i; j++) {
4598 ea[j][i] = ea[i][j];
4604 double fiTerm = TrasanTRK.
fiTerm;
4610 int inlyr(999), outlyr(-1);
4611 int* rStat =
new int[43];
4612 for(
int irStat=0;irStat<43;++irStat)rStat[irStat]=0;
4613 std::vector<MdcRec_wirhit*> pt=TrasanTRK.
hitcol;
4615 if(
debug_ == 4) cout<<
"*********Pt size****"<< pt.size()<<endl;
4618 for (vector<MdcRec_wirhit*>::iterator ii = pt.end()-1;
4619 ii != pt.begin()-1; ii--) {
4620 Num[(*ii)->geo->Lyr()->Id()]++;
4624 for (vector<MdcRec_wirhit*>::iterator ii = pt.end()-1;
4625 ii != pt.begin()-1; ii--) {
4628 if (Num[(*ii)->geo->Lyr()->Id()]>3) {
4630 cout <<
"WARNING: I found " << Num[(*ii)->geo->Lyr()->Id()]
4631 <<
" hits in the layer "
4632 << (*ii)->geo->Lyr()->Id() << std::endl;
4659 double dist[2] = {rechit.
ddl, rechit.
ddr};
4660 double erdist[2] = {rechit.
erddl, rechit.
erddr};
4665 if (rechit.
lr==2 || rechit.
lr==0) lr_decision=-1;
4667 else if (rechit.
lr==1) lr_decision=1;
4670 int ind(geo->
Lyr()->
Id());
4672 lr_decision, rechit.
tdc,
4677 if (inlyr>ind) inlyr = ind;
4678 if (outlyr<ind) outlyr = ind;
4681 cout <<
"**** NUMBER OF Mdc HITS (TRASAN) = " << hit_asso << std::endl;
4684 int empty_between(0), empty(0);
4685 for (
int i= inlyr; i <= outlyr; i++)
4686 if (!rStat[i]) empty_between++;
4687 empty = empty_between+inlyr+(42-outlyr);
4692 track_lead.
type(type);
4693 unsigned int nhit = track_lead.
HitsMdc().size();
4694 if (!nhit &&
debug_ == 4) {
4695 cout <<
" ATTENTION TRACK WITH ONLY HITS " << nhit << std::endl;
4700 double KalFitst(0), KalFitax(0), KalFitschi2(0);
4702 Hep3Vector outer_pivot(track_lead.
x(fiTerm));
4705 std::cout<<
"before track_lead.pivot(outer_pivot) ,the error matrix of track_lead is .."<<track_lead.
Ea()<<std::endl;
4707 track_lead.
pivot(outer_pivot);
4712 HepSymMatrix Eakal(5,0);
4716 double costheta = track_lead.
a()[4] / sqrt(1.0 + track_lead.
a()[4]*track_lead.
a()[4]);
4726 std::cout <<
"from Mdc Pattern Recognition: " << std::endl;
4732 std::cout <<
" dr = " << work.
a()[0]
4733 <<
", Er_dr = " << sqrt(work.
Ea()[0][0]) << std::endl;
4734 std::cout <<
" phi0 = " << work.
a()[1]
4735 <<
", Er_phi0 = " << sqrt(work.
Ea()[1][1]) << std::endl;
4736 std::cout <<
" PT = " << 1/work.
a()[2]
4737 <<
", Er_kappa = " << sqrt(work.
Ea()[2][2]) << std::endl;
4738 std::cout <<
" dz = " << work.
a()[3]
4739 <<
", Er_dz = " << sqrt(work.
Ea()[3][3]) << std::endl;
4740 std::cout <<
" tanl = " << work.
a()[4]
4741 <<
", Er_tanl = " << sqrt(work.
Ea()[4][4]) << std::endl;
4749 cout <<
" Mdc FIRST KALMAN FIT " << std::endl;
4754 cout <<
" dr = " << work.
a()[0]
4755 <<
", Er_dr = " << sqrt(work.
Ea()[0][0]) << std::endl;
4756 cout <<
" phi0 = " << work.
a()[1]
4757 <<
", Er_phi0 = " << sqrt(work.
Ea()[1][1]) << std::endl;
4758 cout <<
" PT = " << 1/work.
a()[2]
4759 <<
", Er_kappa = " << sqrt(work.
Ea()[2][2]) << std::endl;
4760 cout <<
" dz = " << work.
a()[3]
4761 <<
", Er_dz = " << sqrt(work.
Ea()[3][3]) << std::endl;
4762 cout <<
" tanl = " << work.
a()[4]
4763 <<
", Er_tanl = " << sqrt(work.
Ea()[4][4]) << std::endl;
4770 complete_track(TrasanTRK, TrasanTRK_add, track_lead, kaltrk,kalcol,segcol);
4776 DataObject *aRecKalEvent;
4777 eventSvc()->findObject(
"/Event/Recon/RecMdcKalTrackCol", aRecKalEvent);
4778 if(aRecKalEvent!=NULL) {
4780 kalsc = eventSvc()->unregisterObject(
"/Event/Recon/RecMdcKalTrackCol");
4781 if(kalsc != StatusCode::SUCCESS) {
4782 log << MSG::FATAL <<
"Could not unregister RecMdcKalTrack collection" << endreq;
4787 kalsc = eventSvc()->registerObject(
"/Event/Recon/RecMdcKalTrackCol", kalcol);
4788 if( kalsc.isFailure()) {
4789 log << MSG::FATAL <<
"Could not register RecMdcKalTrack" << endreq;
4792 log << MSG::INFO <<
"RecMdcKalTrackCol registered successfully!" <<endreq;
4798 DataObject *aRecKalSegEvent;
4799 eventSvc()->findObject(
"/Event/Recon/RecMdcKalHelixSegCol", aRecKalSegEvent);
4800 if(aRecKalSegEvent!=NULL) {
4802 segsc = eventSvc()->unregisterObject(
"/Event/Recon/RecMdcKalHelixSegCol");
4803 if(segsc != StatusCode::SUCCESS) {
4804 log << MSG::FATAL <<
"Could not unregister RecMdcKalHelixSegCol collection" << endreq;
4809 segsc = eventSvc()->registerObject(
"/Event/Recon/RecMdcKalHelixSegCol", segcol);
4810 if( segsc.isFailure() ) {
4811 log << MSG::FATAL <<
"Could not register RecMdcKalHelixSeg" << endreq;
4814 log << MSG::INFO <<
"RecMdcKalHelixSegCol registered successfully!" <<endreq;
4817 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.);
4818 double r1(0.),r2(0.),kap1(999.),kap2(999.),tanl1(0.),tanl2(0.);
4821 SmartDataPtr<RecMdcKalTrackCol> kaltrkCol(eventSvc(),
"/Event/Recon/RecMdcKalTrackCol");
4823 log << MSG::FATAL <<
"Could not find RecMdcKalTrackCol" << endreq;
4826 log << MSG::INFO <<
"Begin to check RecMdcKalTrackCol"<<endreq;
4827 RecMdcKalTrackCol::iterator iter_trk = kaltrkCol->begin();
4828 for(
int jj=1; iter_trk != kaltrkCol->end(); iter_trk++,jj++) {
4829 log << MSG::DEBUG <<
"retrieved MDC Kalmantrack:"
4830 <<
"Track Id: " << (*iter_trk)->getTrackId()
4831 <<
" Mass of the fit: "<< (*iter_trk)->getMass(2)<< endreq
4832 <<
" Length of the track: "<< (*iter_trk)->getLength(2)
4833 <<
" Tof of the track: "<< (*iter_trk)->getTof(2) << endreq
4834 <<
" Chisq of the fit: "<< (*iter_trk)->getChisq(0,2)
4835 <<
" "<< (*iter_trk)->getChisq(1,2) << endreq
4836 <<
"Ndf of the fit: "<< (*iter_trk)->getNdf(0,1)
4837 <<
" "<< (*iter_trk)->getNdf(1,2) << endreq
4838 <<
"Kappa " << (*iter_trk)->getZHelix()[2]
4843 std::cout<<
"the size of gothelixsegs ..."<<gothelixsegs.size()<<std::endl;
4846 HelixSegRefVec::iterator it_gothelixseg = gothelixsegs.begin();
4847 for( ; it_gothelixseg != gothelixsegs.end(); it_gothelixseg++) {
4849 std::cout<<
"the layerId of this helixseg is ..."<<(*it_gothelixseg)->getLayerId()<<std::endl;
4850 std::cout<<
"the residual of this helixseg exclude the meas hit"<<(*it_gothelixseg)->getResExcl()<<std::endl;
4851 std::cout<<
"the residual of this helixseg include the meas hit"<<(*it_gothelixseg)->getResIncl()<<std::endl;
4852 std::cout<<
"the track id of the helixseg is ..."<<(*it_gothelixseg)->getTrackId() <<std::endl;
4853 std::cout<<
"the tof of the helixseg is ..."<<(*it_gothelixseg)->getTof()<<std::endl;
4854 std::cout<<
"the Zhit of the helixseg is ..."<<(*it_gothelixseg)->getZhit()<<std::endl;
4857 for(
int i = 0; i<43; i++) {
4858 log << MSG::DEBUG <<
"retrieved pathl["<<i<<
"]= "
4859 << (*iter_trk)->getPathl(i) <<endreq;
4863 m_trackid = (*iter_trk)->getTrackId();
4864 for(
int jj =0, iii=0; jj<5; jj++) {
4865 m_length[jj] = (*iter_trk)->getLength(jj);
4866 m_tof[jj] = (*iter_trk)->getTof(jj);
4867 m_nhits[jj] = (*iter_trk)->getNhits(jj);
4868 m_zhelix[jj] = (*iter_trk)->getZHelix()[jj];
4869 m_zhelixe[jj] = (*iter_trk)->getZHelixE()[jj];
4870 m_zhelixmu[jj] = (*iter_trk)->getZHelixMu()[jj];
4871 m_zhelixk[jj] = (*iter_trk)->getZHelixK()[jj];
4872 m_zhelixp[jj] = (*iter_trk)->getZHelixP()[jj];
4873 m_fhelix[jj] = (*iter_trk)->getFHelix()[jj];
4874 m_fhelixe[jj] = (*iter_trk)->getFHelixE()[jj];
4875 m_fhelixmu[jj] = (*iter_trk)->getFHelixMu()[jj];
4876 m_fhelixk[jj] = (*iter_trk)->getFHelixK()[jj];
4877 m_fhelixp[jj] = (*iter_trk)->getFHelixP()[jj];
4878 m_lhelix[jj] = (*iter_trk)->getLHelix()[jj];
4879 m_lhelixe[jj] = (*iter_trk)->getLHelixE()[jj];
4880 m_lhelixmu[jj] = (*iter_trk)->getLHelixMu()[jj];
4881 m_lhelixk[jj] = (*iter_trk)->getLHelixK()[jj];
4882 m_lhelixp[jj] = (*iter_trk)->getLHelixP()[jj];
4884 for(
int kk=0; kk<=jj; kk++,iii++) {
4885 m_zerror[iii] = (*iter_trk)->getZError()[jj][kk];
4886 m_zerrore[iii] = (*iter_trk)->getZErrorE()[jj][kk];
4887 m_zerrormu[iii] = (*iter_trk)->getZErrorMu()[jj][kk];
4888 m_zerrork[iii] = (*iter_trk)->getZErrorK()[jj][kk];
4889 m_zerrorp[iii] = (*iter_trk)->getZErrorP()[jj][kk];
4890 m_ferror[iii] = (*iter_trk)->getFError()[jj][kk];
4891 m_ferrore[iii] = (*iter_trk)->getFErrorE()[jj][kk];
4892 m_ferrormu[iii] = (*iter_trk)->getFErrorMu()[jj][kk];
4893 m_ferrork[iii] = (*iter_trk)->getFErrorK()[jj][kk];
4894 m_ferrorp[iii] = (*iter_trk)->getFErrorP()[jj][kk];
4895 m_lerror[iii] = (*iter_trk)->getLError()[jj][kk];
4896 m_lerrore[iii] = (*iter_trk)->getLErrorE()[jj][kk];
4897 m_lerrormu[iii] = (*iter_trk)->getLErrorMu()[jj][kk];
4898 m_lerrork[iii] = (*iter_trk)->getLErrorK()[jj][kk];
4899 m_lerrorp[iii] = (*iter_trk)->getLErrorP()[jj][kk];
4940 m_chisq[0][0] = (*iter_trk)->getChisq(0,0);
4941 m_chisq[1][0] = (*iter_trk)->getChisq(0,1);
4942 m_chisq[2][0] = (*iter_trk)->getChisq(0,2);
4943 m_chisq[3][0] = (*iter_trk)->getChisq(0,3);
4944 m_chisq[4][0] = (*iter_trk)->getChisq(0,4);
4945 m_chisq[0][1] = (*iter_trk)->getChisq(1,0);
4946 m_chisq[1][1] = (*iter_trk)->getChisq(1,1);
4947 m_chisq[2][1] = (*iter_trk)->getChisq(1,2);
4948 m_chisq[3][1] = (*iter_trk)->getChisq(1,3);
4949 m_chisq[4][1] = (*iter_trk)->getChisq(1,4);
4951 m_ndf[0][0] = (*iter_trk)->getNdf(0,0);
4952 m_ndf[1][0] = (*iter_trk)->getNdf(0,1);
4953 m_ndf[2][0] = (*iter_trk)->getNdf(0,2);
4954 m_ndf[3][0] = (*iter_trk)->getNdf(0,3);
4955 m_ndf[4][0] = (*iter_trk)->getNdf(0,4);
4956 m_ndf[0][1] = (*iter_trk)->getNdf(1,0);
4957 m_ndf[1][1] = (*iter_trk)->getNdf(1,1);
4958 m_ndf[2][1] = (*iter_trk)->getNdf(1,2);
4959 m_ndf[3][1] = (*iter_trk)->getNdf(1,3);
4960 m_ndf[4][1] = (*iter_trk)->getNdf(1,4);
4962 m_stat[0][0] = (*iter_trk)->getStat(0,0);
4963 m_stat[1][0] = (*iter_trk)->getStat(0,1);
4964 m_stat[2][0] = (*iter_trk)->getStat(0,2);
4965 m_stat[3][0] = (*iter_trk)->getStat(0,3);
4966 m_stat[4][0] = (*iter_trk)->getStat(0,4);
4967 m_stat[0][1] = (*iter_trk)->getStat(1,0);
4968 m_stat[1][1] = (*iter_trk)->getStat(1,1);
4969 m_stat[2][1] = (*iter_trk)->getStat(1,2);
4970 m_stat[3][1] = (*iter_trk)->getStat(1,3);
4971 m_stat[4][1] = (*iter_trk)->getStat(1,4);
4973 m_fptot = sqrt(1+pow(m_fhelix[4],2))/m_fhelix[2];
4974 m_fptote = sqrt(1+pow(m_fhelixe[4],2))/m_fhelixe[2];
4975 m_fptotmu = sqrt(1+pow(m_fhelixmu[4],2))/m_fhelixmu[2];
4976 m_fptotk = sqrt(1+pow(m_fhelixk[4],2))/m_fhelixk[2];
4977 m_fptotp = sqrt(1+pow(m_fhelixp[4],2))/m_fhelixp[2];
4979 m_zpt = 1/m_zhelix[2];
4980 m_zpte = 1/m_zhelixe[2];
4981 m_zptmu = 1/m_zhelixmu[2];
4982 m_zptk = 1/m_zhelixk[2];
4983 m_zptp = 1/m_zhelixp[2];
4985 m_fpt = 1/m_fhelix[2];
4986 m_fpte = 1/m_fhelixe[2];
4987 m_fptmu = 1/m_fhelixmu[2];
4988 m_fptk = 1/m_fhelixk[2];
4989 m_fptp = 1/m_fhelixp[2];
4991 m_lpt = 1/m_lhelix[2];
4992 m_lpte = 1/m_lhelixe[2];
4993 m_lptmu = 1/m_lhelixmu[2];
4994 m_lptk = 1/m_lhelixk[2];
4995 m_lptp = 1/m_lhelixp[2];
4997 m_lptot = sqrt(1+pow(m_lhelix[4],2))/m_lhelix[2];
4998 m_lptote = sqrt(1+pow(m_lhelixe[4],2))/m_lhelixe[2];
4999 m_lptotmu = sqrt(1+pow(m_lhelixmu[4],2))/m_lhelixmu[2];
5000 m_lptotk = sqrt(1+pow(m_lhelixk[4],2))/m_lhelixk[2];
5001 m_lptotp = sqrt(1+pow(m_lhelixp[4],2))/m_lhelixp[2];
5003 m_zptot = sqrt(1+pow(m_zhelix[4],2))/m_zhelix[2];
5004 m_zptote = sqrt(1+pow(m_zhelixe[4],2))/m_zhelixe[2];
5005 m_zptotmu = sqrt(1+pow(m_zhelixmu[4],2))/m_zhelixmu[2];
5006 m_zptotk = sqrt(1+pow(m_zhelixk[4],2))/m_zhelixk[2];
5007 m_zptotp = sqrt(1+pow(m_zhelixp[4],2))/m_zhelixp[2];
5009 m_zsigp = sqrt(pow((m_zptot/m_zhelix[2]),2)*m_zerror[5]+
5010 pow((m_zhelix[4]/m_zptot),2)*pow((1/m_zhelix[2]),4)*m_zerror[14]-
5011 2*m_zhelix[4]*m_zerror[12]*pow((1/m_zhelix[2]),3));
5012 m_zsigpe = sqrt(pow((m_zptote/m_zhelixe[2]),2)*m_zerrore[5]+
5013 pow((m_zhelixe[4]/m_zptote),2)*pow((1/m_zhelixe[2]),4)*m_zerrore[14]-
5014 2*m_zhelixe[4]*m_zerrore[12]*pow((1/m_zhelixe[2]),3));
5015 m_zsigpmu = sqrt(pow((m_zptotmu/m_zhelixmu[2]),2)*m_zerrormu[5]+
5016 pow((m_zhelixmu[4]/m_zptotmu),2)*pow((1/m_zhelixmu[2]),4)*m_zerrormu[14]-
5017 2*m_zhelixmu[4]*m_zerrormu[12]*pow((1/m_zhelixmu[2]),3));
5018 m_zsigpk = sqrt(pow((m_zptotk/m_zhelixk[2]),2)*m_zerrork[5]+
5019 pow((m_zhelixk[4]/m_zptotk),2)*pow((1/m_zhelixk[2]),4)*m_zerrork[14]-
5020 2*m_zhelixk[4]*m_zerrork[12]*pow((1/m_zhelixk[2]),3));
5021 m_zsigpp = sqrt(pow((m_zptotp/m_zhelixp[2]),2)*m_zerrorp[5]+
5022 pow((m_zhelixp[4]/m_zptotp),2)*pow((1/m_zhelixp[2]),4)*m_zerrorp[14]-
5023 2*m_zhelixp[4]*m_zerrorp[12]*pow((1/m_zhelixp[2]),3));
5026 StatusCode sc1 = m_nt1->write();
5027 if( sc1.isFailure() ) cout<<
"Ntuple1 filling failed!"<<endl;
5032 phi1 = (*iter_trk)->getFFi0();
5033 r1 = (*iter_trk)->getFDr();
5034 z1 = (*iter_trk)->getFDz();
5035 kap1 = (*iter_trk)->getFCpa();
5036 tanl1 = (*iter_trk)->getFTanl();
5039 p1 = sqrt(1+tanl1*tanl1)/kap1;
5040 the1 =
M_PI/2-atan(tanl1);
5041 }
else if(jj == 2) {
5042 phi2 = (*iter_trk)->getFFi0();
5043 r2 = (*iter_trk)->getFDr();
5044 z2 = (*iter_trk)->getFDz();
5045 kap2 = (*iter_trk)->getFCpa();
5046 tanl2 = (*iter_trk)->getFTanl();
5049 p2 = sqrt(1+tanl2*tanl2)/kap1;
5050 the2 =
M_PI/2-atan(tanl2);
5058 m_delthe = the1 + the2;
5061 StatusCode sc2 = m_nt2->write();
5062 if( sc2.isFailure() ) cout<<
"Ntuple2 filling failed!"<<endl;
5070 cout <<
"Kalfitting finished " << std::endl;
5077 MsgStream log(
msgSvc(), name());
5078 double Pt_threshold(0.3);
5079 Hep3Vector IP(0,0,0);
5086 if ( !&whMgr )
return;
5089 int ntrk = mdcMgr->size();
5092 int nhits = whMgr->size();
5097 DataObject *aReconEvent;
5098 eventSvc()->findObject(
"/Event/Recon",aReconEvent);
5102 StatusCode sc = eventSvc()->registerObject(
"/Event/Recon",recevt );
5103 if(sc!=StatusCode::SUCCESS) {
5104 log << MSG::FATAL <<
"Could not register ReconEvent" <<endreq;
5112 log << MSG::INFO <<
"beginning to make RecMdcKalTrackCol" <<endreq;
5124 if ((TrasanTRK_add.
decision & 32) == 32 ||
5125 (TrasanTRK_add.
decision & 64) == 64) type = 1;
5130 TrasanTRK.
pivot[2]);
5133 for(
int i = 0; i < 5; i++)
5134 a[i] = TrasanTRK.
helix[i];
5137 for(
int i = 0, k = 0; i < 5; i++) {
5138 for(
int j = 0; j <= i; j++) {
5140 ea[j][i] = ea[i][j];
5146 double fiTerm = TrasanTRK.
fiTerm;
5154 int trasqual = TrasanTRK_add.
quality;
5155 if(
debug_ == 4) cout<<
"kalman_fitting>trasqual... "<<trasqual<<endl;
5156 if (trasqual)
return;
5158 int inlyr(999), outlyr(-1);
5159 int* rStat =
new int[43];
5160 for(
int irStat=0;irStat<43;++irStat)rStat[irStat]=0;
5161 std::vector<MdcRec_wirhit*> pt=TrasanTRK.
hitcol;
5163 if(
debug_ == 4) cout<<
"*********Pt size****"<< pt.size()<<endl;
5166 for (vector<MdcRec_wirhit*>::iterator ii = pt.end()-1;
5167 ii != pt.begin()-1; ii--) {
5168 Num[(*ii)->geo->Lyr()->Id()]++;
5171 for (vector<MdcRec_wirhit*>::iterator ii = pt.end()-1;
5172 ii != pt.begin()-1; ii--) {
5175 if (Num[(*ii)->geo->Lyr()->Id()]>3) {
5177 cout <<
"WARNING: I found " << Num[(*ii)->geo->Lyr()->Id()]
5178 <<
" hits in the layer "
5179 << (*ii)->geo->Lyr()->Id() << std::endl;
5185 double dist[2] = {rechit.
ddl, rechit.
ddr};
5186 double erdist[2] = {rechit.
erddl, rechit.
erddr};
5191 if (rechit.
lr==2 || rechit.
lr==0) lr_decision=-1;
5193 else if (rechit.
lr==1) lr_decision=1;
5196 int ind(geo->
Lyr()->
Id());
5198 lr_decision, rechit.
tdc,
5203 if (inlyr>ind) inlyr = ind;
5204 if (outlyr<ind) outlyr = ind;
5207 int empty_between(0), empty(0);
5208 for (
int i= inlyr; i <= outlyr; i++)
5209 if (!rStat[i]) empty_between++;
5210 empty = empty_between+inlyr+(42-outlyr);
5214 cout <<
"**** NUMBER OF Mdc HITS (TRASAN) = " << hit_asso << std::endl;
5219 track_lead.
type(type);
5220 unsigned int nhit = track_lead.
HitsMdc().size();
5222 cout <<
" ATTENTION TRACK WITH ONLY HITS " << nhit << std::endl;
5227 double KalFitst(0), KalFitax(0), KalFitschi2(0);
5229 Hep3Vector outer_pivot(track_lead.
x(fiTerm));
5232 std::cout<<
"before track_lead.pivot(outer_pivot) ,the error matrix of track_lead is .."<<track_lead.
Ea()<<std::endl;
5234 track_lead.
pivot(outer_pivot);
5239 HepSymMatrix Eakal(5,0);
5243 double costheta = track_lead.
a()[4] / sqrt(1.0 + track_lead.
a()[4]*track_lead.
a()[4]);
5253 std::cout <<
"from Mdc Pattern Recognition: " << std::endl;
5259 std::cout <<
" dr = " << work.
a()[0]
5260 <<
", Er_dr = " << sqrt(work.
Ea()[0][0]) << std::endl;
5261 std::cout <<
" phi0 = " << work.
a()[1]
5262 <<
", Er_phi0 = " << sqrt(work.
Ea()[1][1]) << std::endl;
5263 std::cout <<
" PT = " << 1/work.
a()[2]
5264 <<
", Er_kappa = " << sqrt(work.
Ea()[2][2]) << std::endl;
5265 std::cout <<
" dz = " << work.
a()[3]
5266 <<
", Er_dz = " << sqrt(work.
Ea()[3][3]) << std::endl;
5267 std::cout <<
" tanl = " << work.
a()[4]
5268 <<
", Er_tanl = " << sqrt(work.
Ea()[4][4]) << std::endl;
5275 cout <<
" Mdc FIRST KALMAN FIT " << std::endl;
5280 cout <<
" dr = " << work1.
a()[0]
5281 <<
", Er_dr = " << sqrt(work1.
Ea()[0][0]) << std::endl;
5282 cout <<
" phi0 = " << work1.
a()[1]
5283 <<
", Er_phi0 = " << sqrt(work1.
Ea()[1][1]) << std::endl;
5284 cout <<
" PT = " << 1/work1.
a()[2]
5285 <<
", Er_kappa = " << sqrt(work1.
Ea()[2][2]) << std::endl;
5286 cout <<
" dz = " << work1.
a()[3]
5287 <<
", Er_dz = " << sqrt(work1.
Ea()[3][3]) << std::endl;
5288 cout <<
" tanl = " << work1.
a()[4]
5289 <<
", Er_tanl = " << sqrt(work1.
Ea()[4][4]) << std::endl;
5296 complete_track(TrasanTRK, TrasanTRK_add, track_lead, kaltrk,kalcol,segcol);
5301 DataObject *aRecKalEvent;
5302 eventSvc()->findObject(
"/Event/Recon/RecMdcKalTrackCol", aRecKalEvent);
5303 if(aRecKalEvent!=NULL) {
5305 kalsc = eventSvc()->unregisterObject(
"/Event/Recon/RecMdcKalTrackCol");
5306 if(kalsc != StatusCode::SUCCESS) {
5307 log << MSG::FATAL <<
"Could not unregister RecMdcKalTrack collection" << endreq;
5312 kalsc = eventSvc()->registerObject(
"/Event/Recon/RecMdcKalTrackCol", kalcol);
5313 if( kalsc.isFailure()) {
5314 log << MSG::FATAL <<
"Could not register RecMdcKalTrack" << endreq;
5317 log << MSG::INFO <<
"RecMdcKalTrackCol registered successfully!" <<endreq;
5323 DataObject *aRecKalSegEvent;
5324 eventSvc()->findObject(
"/Event/Recon/RecMdcKalHelixSegCol", aRecKalSegEvent);
5325 if(aRecKalSegEvent!=NULL) {
5327 segsc = eventSvc()->unregisterObject(
"/Event/Recon/RecMdcKalHelixSegCol");
5328 if(segsc != StatusCode::SUCCESS) {
5329 log << MSG::FATAL <<
"Could not unregister RecMdcKalHelixSegCol collection" << endreq;
5334 segsc = eventSvc()->registerObject(
"/Event/Recon/RecMdcKalHelixSegCol", segcol);
5335 if( segsc.isFailure() ) {
5336 log << MSG::FATAL <<
"Could not register RecMdcKalHelixSeg" << endreq;
5339 log << MSG::INFO <<
"RecMdcKalHelixSegCol registered successfully!" <<endreq;
5342 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.);
5343 double r1(0.),r2(0.),kap1(999.),kap2(999.),tanl1(0.),tanl2(0.);
5346 SmartDataPtr<RecMdcKalTrackCol> kaltrkCol(eventSvc(),
"/Event/Recon/RecMdcKalTrackCol");
5348 log << MSG::FATAL <<
"Could not find RecMdcKalTrackCol" << endreq;
5351 log << MSG::INFO <<
"Begin to check RecMdcKalTrackCol"<<endreq;
5352 RecMdcKalTrackCol::iterator iter_trk = kaltrkCol->begin();
5353 for(
int jj=1; iter_trk != kaltrkCol->end(); iter_trk++,jj++) {
5354 log << MSG::DEBUG <<
"retrieved MDC Kalmantrack:"
5355 <<
"Track Id: " << (*iter_trk)->getTrackId()
5356 <<
" Mass of the fit: "<< (*iter_trk)->getMass(2)<< endreq
5357 <<
" Length of the track: "<< (*iter_trk)->getLength(2)
5358 <<
" Tof of the track: "<< (*iter_trk)->getTof(2) << endreq
5359 <<
" Chisq of the fit: "<< (*iter_trk)->getChisq(0,2)
5360 <<
" "<< (*iter_trk)->getChisq(1,2) << endreq
5361 <<
"Ndf of the fit: "<< (*iter_trk)->getNdf(0,1)
5362 <<
" "<< (*iter_trk)->getNdf(1,2) << endreq
5363 <<
"Kappa " << (*iter_trk)->getZHelix()[2]
5364 <<
"zhelixmu "<<(*iter_trk)->getZHelixMu()
5369 std::cout<<
"the size of gothelixsegs ..."<<gothelixsegs.size()<<std::endl;
5372 HelixSegRefVec::iterator it_gothelixseg = gothelixsegs.begin();
5373 for( ; it_gothelixseg != gothelixsegs.end(); it_gothelixseg++) {
5375 std::cout<<
"the layerId of this helixseg is ..."<<(*it_gothelixseg)->getLayerId()<<std::endl;
5376 std::cout<<
"the residual of this helixseg exclude the meas hit"<<(*it_gothelixseg)->getResExcl()<<std::endl;
5377 std::cout<<
"the residual of this helixseg include the meas hit"<<(*it_gothelixseg)->getResIncl()<<std::endl;
5378 std::cout<<
"the track id of the helixseg is ..."<<(*it_gothelixseg)->getTrackId() <<std::endl;
5379 std::cout<<
"the tof of the helixseg is ..."<<(*it_gothelixseg)->getTof()<<std::endl;
5380 std::cout<<
"the Zhit of the helixseg is ..."<<(*it_gothelixseg)->getZhit()<<std::endl;
5383 for(
int i = 0; i<43; i++) {
5384 log << MSG::DEBUG <<
"retrieved pathl["<<i<<
"]= "
5385 << (*iter_trk)->getPathl(i) <<endreq;
5389 m_trackid = (*iter_trk)->getTrackId();
5390 for(
int jj =0, iii=0; jj<5; jj++) {
5391 m_length[jj] = (*iter_trk)->getLength(jj);
5392 m_tof[jj] = (*iter_trk)->getTof(jj);
5393 m_nhits[jj] = (*iter_trk)->getNhits(jj);
5394 m_zhelix[jj] = (*iter_trk)->getZHelix()[jj];
5395 m_zhelixe[jj] = (*iter_trk)->getZHelixE()[jj];
5396 m_zhelixmu[jj] = (*iter_trk)->getZHelixMu()[jj];
5397 m_zhelixk[jj] = (*iter_trk)->getZHelixK()[jj];
5398 m_zhelixp[jj] = (*iter_trk)->getZHelixP()[jj];
5399 m_fhelix[jj] = (*iter_trk)->getFHelix()[jj];
5400 m_fhelixe[jj] = (*iter_trk)->getFHelixE()[jj];
5401 m_fhelixmu[jj] = (*iter_trk)->getFHelixMu()[jj];
5402 m_fhelixk[jj] = (*iter_trk)->getFHelixK()[jj];
5403 m_fhelixp[jj] = (*iter_trk)->getFHelixP()[jj];
5404 m_lhelix[jj] = (*iter_trk)->getLHelix()[jj];
5405 m_lhelixe[jj] = (*iter_trk)->getLHelixE()[jj];
5406 m_lhelixmu[jj] = (*iter_trk)->getLHelixMu()[jj];
5407 m_lhelixk[jj] = (*iter_trk)->getLHelixK()[jj];
5408 m_lhelixp[jj] = (*iter_trk)->getLHelixP()[jj];
5410 for(
int kk=0; kk<=jj; kk++,iii++) {
5411 m_zerror[iii] = (*iter_trk)->getZError()[jj][kk];
5412 m_zerrore[iii] = (*iter_trk)->getZErrorE()[jj][kk];
5413 m_zerrormu[iii] = (*iter_trk)->getZErrorMu()[jj][kk];
5414 m_zerrork[iii] = (*iter_trk)->getZErrorK()[jj][kk];
5415 m_zerrorp[iii] = (*iter_trk)->getZErrorP()[jj][kk];
5416 m_ferror[iii] = (*iter_trk)->getFError()[jj][kk];
5417 m_ferrore[iii] = (*iter_trk)->getFErrorE()[jj][kk];
5418 m_ferrormu[iii] = (*iter_trk)->getFErrorMu()[jj][kk];
5419 m_ferrork[iii] = (*iter_trk)->getFErrorK()[jj][kk];
5420 m_ferrorp[iii] = (*iter_trk)->getFErrorP()[jj][kk];
5421 m_lerror[iii] = (*iter_trk)->getLError()[jj][kk];
5422 m_lerrore[iii] = (*iter_trk)->getLErrorE()[jj][kk];
5423 m_lerrormu[iii] = (*iter_trk)->getLErrorMu()[jj][kk];
5424 m_lerrork[iii] = (*iter_trk)->getLErrorK()[jj][kk];
5425 m_lerrorp[iii] = (*iter_trk)->getLErrorP()[jj][kk];
5432 m_chisq[0][0] = (*iter_trk)->getChisq(0,0);
5433 m_chisq[1][0] = (*iter_trk)->getChisq(0,1);
5434 m_chisq[2][0] = (*iter_trk)->getChisq(0,2);
5435 m_chisq[3][0] = (*iter_trk)->getChisq(0,3);
5436 m_chisq[4][0] = (*iter_trk)->getChisq(0,4);
5437 m_chisq[0][1] = (*iter_trk)->getChisq(1,0);
5438 m_chisq[1][1] = (*iter_trk)->getChisq(1,1);
5439 m_chisq[2][1] = (*iter_trk)->getChisq(1,2);
5440 m_chisq[3][1] = (*iter_trk)->getChisq(1,3);
5441 m_chisq[4][1] = (*iter_trk)->getChisq(1,4);
5443 m_ndf[0][0] = (*iter_trk)->getNdf(0,0);
5444 m_ndf[1][0] = (*iter_trk)->getNdf(0,1);
5445 m_ndf[2][0] = (*iter_trk)->getNdf(0,2);
5446 m_ndf[3][0] = (*iter_trk)->getNdf(0,3);
5447 m_ndf[4][0] = (*iter_trk)->getNdf(0,4);
5448 m_ndf[0][1] = (*iter_trk)->getNdf(1,0);
5449 m_ndf[1][1] = (*iter_trk)->getNdf(1,1);
5450 m_ndf[2][1] = (*iter_trk)->getNdf(1,2);
5451 m_ndf[3][1] = (*iter_trk)->getNdf(1,3);
5452 m_ndf[4][1] = (*iter_trk)->getNdf(1,4);
5454 m_stat[0][0] = (*iter_trk)->getStat(0,0);
5455 m_stat[1][0] = (*iter_trk)->getStat(0,1);
5456 m_stat[2][0] = (*iter_trk)->getStat(0,2);
5457 m_stat[3][0] = (*iter_trk)->getStat(0,3);
5458 m_stat[4][0] = (*iter_trk)->getStat(0,4);
5459 m_stat[0][1] = (*iter_trk)->getStat(1,0);
5460 m_stat[1][1] = (*iter_trk)->getStat(1,1);
5461 m_stat[2][1] = (*iter_trk)->getStat(1,2);
5462 m_stat[3][1] = (*iter_trk)->getStat(1,3);
5463 m_stat[4][1] = (*iter_trk)->getStat(1,4);
5465 m_fptot = sqrt(1+pow(m_fhelix[4],2))/m_fhelix[2];
5466 m_fptote = sqrt(1+pow(m_fhelixe[4],2))/m_fhelixe[2];
5467 m_fptotmu = sqrt(1+pow(m_fhelixmu[4],2))/m_fhelixmu[2];
5468 m_fptotk = sqrt(1+pow(m_fhelixk[4],2))/m_fhelixk[2];
5469 m_fptotp = sqrt(1+pow(m_fhelixp[4],2))/m_fhelixp[2];
5471 m_zpt = 1/m_zhelix[2];
5472 m_zpte = 1/m_zhelixe[2];
5473 m_zptmu = 1/m_zhelixmu[2];
5474 m_zptk = 1/m_zhelixk[2];
5475 m_zptp = 1/m_zhelixp[2];
5477 m_fpt = 1/m_fhelix[2];
5478 m_fpte = 1/m_fhelixe[2];
5479 m_fptmu = 1/m_fhelixmu[2];
5480 m_fptk = 1/m_fhelixk[2];
5481 m_fptp = 1/m_fhelixp[2];
5483 m_lpt = 1/m_lhelix[2];
5484 m_lpte = 1/m_lhelixe[2];
5485 m_lptmu = 1/m_lhelixmu[2];
5486 m_lptk = 1/m_lhelixk[2];
5487 m_lptp = 1/m_lhelixp[2];
5489 m_lptot = sqrt(1+pow(m_lhelix[4],2))/m_lhelix[2];
5490 m_lptote = sqrt(1+pow(m_lhelixe[4],2))/m_lhelixe[2];
5491 m_lptotmu = sqrt(1+pow(m_lhelixmu[4],2))/m_lhelixmu[2];
5492 m_lptotk = sqrt(1+pow(m_lhelixk[4],2))/m_lhelixk[2];
5493 m_lptotp = sqrt(1+pow(m_lhelixp[4],2))/m_lhelixp[2];
5495 m_zptot = sqrt(1+pow(m_zhelix[4],2))/m_zhelix[2];
5496 m_zptote = sqrt(1+pow(m_zhelixe[4],2))/m_zhelixe[2];
5497 m_zptotmu = sqrt(1+pow(m_zhelixmu[4],2))/m_zhelixmu[2];
5498 m_zptotk = sqrt(1+pow(m_zhelixk[4],2))/m_zhelixk[2];
5499 m_zptotp = sqrt(1+pow(m_zhelixp[4],2))/m_zhelixp[2];
5501 m_zsigp = sqrt(pow((m_zptot/m_zhelix[2]),2)*m_zerror[5]+
5502 pow((m_zhelix[4]/m_zptot),2)*pow((1/m_zhelix[2]),4)*m_zerror[14]-
5503 2*m_zhelix[4]*m_zerror[12]*pow((1/m_zhelix[2]),3));
5504 m_zsigpe = sqrt(pow((m_zptote/m_zhelixe[2]),2)*m_zerrore[5]+
5505 pow((m_zhelixe[4]/m_zptote),2)*pow((1/m_zhelixe[2]),4)*m_zerrore[14]-
5506 2*m_zhelixe[4]*m_zerrore[12]*pow((1/m_zhelixe[2]),3));
5507 m_zsigpmu = sqrt(pow((m_zptotmu/m_zhelixmu[2]),2)*m_zerrormu[5]+
5508 pow((m_zhelixmu[4]/m_zptotmu),2)*pow((1/m_zhelixmu[2]),4)*m_zerrormu[14]-
5509 2*m_zhelixmu[4]*m_zerrormu[12]*pow((1/m_zhelixmu[2]),3));
5510 m_zsigpk = sqrt(pow((m_zptotk/m_zhelixk[2]),2)*m_zerrork[5]+
5511 pow((m_zhelixk[4]/m_zptotk),2)*pow((1/m_zhelixk[2]),4)*m_zerrork[14]-
5512 2*m_zhelixk[4]*m_zerrork[12]*pow((1/m_zhelixk[2]),3));
5513 m_zsigpp = sqrt(pow((m_zptotp/m_zhelixp[2]),2)*m_zerrorp[5]+
5514 pow((m_zhelixp[4]/m_zptotp),2)*pow((1/m_zhelixp[2]),4)*m_zerrorp[14]-
5515 2*m_zhelixp[4]*m_zerrorp[12]*pow((1/m_zhelixp[2]),3));
5518 StatusCode sc1 = m_nt1->write();
5519 if( sc1.isFailure() ) cout<<
"Ntuple1 filling failed!"<<endl;
5524 phi1 = (*iter_trk)->getFFi0();
5525 r1 = (*iter_trk)->getFDr();
5526 z1 = (*iter_trk)->getFDz();
5527 kap1 = (*iter_trk)->getFCpa();
5528 tanl1 = (*iter_trk)->getFTanl();
5531 p1 = sqrt(1+tanl1*tanl1)/kap1;
5532 the1 =
M_PI/2-atan(tanl1);
5533 }
else if(jj == 2) {
5534 phi2 = (*iter_trk)->getFFi0();
5535 r2 = (*iter_trk)->getFDr();
5536 z2 = (*iter_trk)->getFDz();
5537 kap2 = (*iter_trk)->getFCpa();
5538 tanl2 = (*iter_trk)->getFTanl();
5541 p2 = sqrt(1+tanl2*tanl2)/kap1;
5542 the2 =
M_PI/2-atan(tanl2);
5550 m_delthe = the1 + the2;
5553 StatusCode sc2 = m_nt2->write();
5554 if( sc2.isFailure() ) cout<<
"Ntuple2 filling failed!"<<endl;
5558 cout <<
"Kalfitting finished " << std::endl;
5566 MsgStream log(
msgSvc(), name());
5567 double Pt_threshold(0.3);
5568 Hep3Vector IP(0,0,0);
5575 if ( !&whMgr )
return;
5578 int ntrk = mdcMgr->size();
5581 int nhits = whMgr->size();
5585 double* rY =
new double[ntrk];
5586 double* rfiTerm =
new double[ntrk];
5587 double* rPt =
new double[ntrk];
5588 int* rOM =
new int[ntrk];
5589 unsigned int* order =
new unsigned int[ntrk];
5590 unsigned int* rCont =
new unsigned int[ntrk];
5591 unsigned int* rGen =
new unsigned int[ntrk];
5594 Hep3Vector csmp3[2];
5595 for(vector<MdcRec_trk>::iterator it = mdcMgr->begin(),
5596 end = mdcMgr->end(); it != end; it++) {
5598 rfiTerm[index]=it->fiTerm;
5603 rPt[index] = 1 / fabs(it->helix[2]);
5604 if(
debug_ == 4) cout<<
"rPt...[ "<<index<<
" ]...."<< rPt[index] <<endl;
5605 if(rPt[index] < 0) rPt[index] =
DBL_MAX;
5607 std::vector<MdcRec_wirhit*> pt = it->hitcol ;
5608 if(
debug_ == 4) cout<<
"ppt size: "<< pt.size()<<endl;
5610 for (vector<MdcRec_wirhit*>::iterator ii = pt.end()-1;
5611 ii !=pt.begin()-1; ii--) {
5612 int lyr((*ii)->geo->Lyr()->Id());
5613 if (outermost < lyr) {
5615 rY[index] = (*ii)->geo->Forward().y();
5617 if(
debug_ == 4) cout<<
"outmost: "<<outermost<<
" lyr: "<<lyr<<endl;
5619 rOM[index] = outermost;
5620 order[index] = index;
5625 for (
int j, k = ntrk - 1; k >= 0; k = j){
5627 for(
int i = 1; i <= k; i++)
5628 if(rY[i - 1] < rY[i]){
5630 std::swap(order[i], order[j]);
5631 std::swap(rY[i], rY[j]);
5632 std::swap(rOM[i], rOM[j]);
5633 std::swap(rCont[i], rCont[j]);
5634 std::swap(rGen[i], rGen[j]);
5643 DataObject *aReconEvent;
5644 eventSvc()->findObject(
"/Event/Recon",aReconEvent);
5648 StatusCode sc = eventSvc()->registerObject(
"/Event/Recon",recevt );
5649 if(sc!=StatusCode::SUCCESS) {
5650 log << MSG::FATAL <<
"Could not register ReconEvent" <<endreq;
5658 log << MSG::INFO <<
"beginning to make RecMdcKalTrackCol" <<endreq;
5665 MdcRec_trk& TrasanTRK = *(mdcMgr->begin() + order[1]);
5674 cout <<
"******* KalFit NUMBER : " << newcount << std::endl;
5678 if ((TrasanTRK_add.
decision & 32) == 32 ||
5679 (TrasanTRK_add.
decision & 64) == 64) type = 1;
5684 TrasanTRK.
pivot[2]);
5687 for(
int i = 0; i < 5; i++)
5688 a[i] = TrasanTRK.
helix[i];
5691 for(
int i = 0, k = 0; i < 5; i++) {
5692 for(
int j = 0; j <= i; j++) {
5694 ea[j][i] = ea[i][j];
5700 double fiTerm = TrasanTRK.
fiTerm;
5707 for(
int l = 0; l < ntrk; l++) {
5708 MdcRec_trk& TrasanTRK1 = *(mdcMgr->begin() + order[l]);
5709 MdcRec_trk_add& TrasanTRK_add1 = *(mdc_addMgr->begin()+order[l]);
5711 int trasqual = TrasanTRK_add1.
quality;
5712 if(
debug_ == 4) cout<<
"kalman_fitting>trasqual... "<<trasqual<<endl;
5713 if (trasqual)
continue;
5715 int inlyr(999), outlyr(-1);
5716 int* rStat =
new int[43];
5717 for(
int irStat=0;irStat<43;++irStat)rStat[irStat]=0;
5718 std::vector<MdcRec_wirhit*> pt=TrasanTRK1.
hitcol;
5720 if(
debug_ == 4) cout<<
"*********Pt size****"<< pt.size()<<endl;
5723 for (vector<MdcRec_wirhit*>::iterator ii = pt.end()-1;
5724 ii != pt.begin()-1; ii--) {
5725 Num[(*ii)->geo->Lyr()->Id()]++;
5728 for (vector<MdcRec_wirhit*>::iterator ii = pt.end()-1;
5729 ii != pt.begin()-1; ii--) {
5732 if (Num[(*ii)->geo->Lyr()->Id()]>3) {
5734 cout <<
"WARNING: I found " << Num[(*ii)->geo->Lyr()->Id()]
5735 <<
" hits in the layer "
5736 << (*ii)->geo->Lyr()->Id() << std::endl;
5742 double dist[2] = {rechit.
ddl, rechit.
ddr};
5743 double erdist[2] = {rechit.
erddl, rechit.
erddr};
5748 if (rechit.
lr==2 || rechit.
lr==0) lr_decision=-1;
5750 else if (rechit.
lr==1) lr_decision=1;
5753 int ind(geo->
Lyr()->
Id());
5755 lr_decision, rechit.
tdc,
5760 if (inlyr>ind) inlyr = ind;
5761 if (outlyr<ind) outlyr = ind;
5764 int empty_between(0), empty(0);
5765 for (
int i= inlyr; i <= outlyr; i++)
5766 if (!rStat[i]) empty_between++;
5767 empty = empty_between+inlyr+(42-outlyr);
5771 cout <<
"**** NUMBER OF Mdc HITS (TRASAN) = " << hit_asso << std::endl;
5776 track_lead.
type(type);
5777 unsigned int nhit = track_lead.
HitsMdc().size();
5779 cout <<
" ATTENTION TRACK WITH ONLY HITS " << nhit << std::endl;
5784 double KalFitst(0), KalFitax(0), KalFitschi2(0);
5786 Hep3Vector outer_pivot(track_lead.
x(fiTerm));
5789 std::cout<<
"before track_lead.pivot(outer_pivot) ,the error matrix of track_lead is .."<<track_lead.
Ea()<<std::endl;
5791 track_lead.
pivot(outer_pivot);
5796 HepSymMatrix Eakal(5,0);
5800 double costheta = track_lead.
a()[4] / sqrt(1.0 + track_lead.
a()[4]*track_lead.
a()[4]);
5810 std::cout <<
"from Mdc Pattern Recognition: " << std::endl;
5816 std::cout <<
" dr = " << work.
a()[0]
5817 <<
", Er_dr = " << sqrt(work.
Ea()[0][0]) << std::endl;
5818 std::cout <<
" phi0 = " << work.
a()[1]
5819 <<
", Er_phi0 = " << sqrt(work.
Ea()[1][1]) << std::endl;
5820 std::cout <<
" PT = " << 1/work.
a()[2]
5821 <<
", Er_kappa = " << sqrt(work.
Ea()[2][2]) << std::endl;
5822 std::cout <<
" dz = " << work.
a()[3]
5823 <<
", Er_dz = " << sqrt(work.
Ea()[3][3]) << std::endl;
5824 std::cout <<
" tanl = " << work.
a()[4]
5825 <<
", Er_tanl = " << sqrt(work.
Ea()[4][4]) << std::endl;
5832 cout <<
" Mdc FIRST KALMAN FIT " << std::endl;
5837 cout <<
" dr = " << work1.
a()[0]
5838 <<
", Er_dr = " << sqrt(work1.
Ea()[0][0]) << std::endl;
5839 cout <<
" phi0 = " << work1.
a()[1]
5840 <<
", Er_phi0 = " << sqrt(work1.
Ea()[1][1]) << std::endl;
5841 cout <<
" PT = " << 1/work1.
a()[2]
5842 <<
", Er_kappa = " << sqrt(work1.
Ea()[2][2]) << std::endl;
5843 cout <<
" dz = " << work1.
a()[3]
5844 <<
", Er_dz = " << sqrt(work1.
Ea()[3][3]) << std::endl;
5845 cout <<
" tanl = " << work1.
a()[4]
5846 <<
", Er_tanl = " << sqrt(work1.
Ea()[4][4]) << std::endl;
5853 complete_track(TrasanTRK, TrasanTRK_add, track_lead, kaltrk,kalcol,segcol);
5859 DataObject *aRecKalEvent;
5860 eventSvc()->findObject(
"/Event/Recon/RecMdcKalTrackCol", aRecKalEvent);
5861 if(aRecKalEvent!=NULL) {
5863 kalsc = eventSvc()->unregisterObject(
"/Event/Recon/RecMdcKalTrackCol");
5864 if(kalsc != StatusCode::SUCCESS) {
5865 log << MSG::FATAL <<
"Could not unregister RecMdcKalTrack collection" << endreq;
5870 kalsc = eventSvc()->registerObject(
"/Event/Recon/RecMdcKalTrackCol", kalcol);
5871 if( kalsc.isFailure()) {
5872 log << MSG::FATAL <<
"Could not register RecMdcKalTrack" << endreq;
5875 log << MSG::INFO <<
"RecMdcKalTrackCol registered successfully!" <<endreq;
5881 DataObject *aRecKalSegEvent;
5882 eventSvc()->findObject(
"/Event/Recon/RecMdcKalHelixSegCol", aRecKalSegEvent);
5883 if(aRecKalSegEvent!=NULL) {
5885 segsc = eventSvc()->unregisterObject(
"/Event/Recon/RecMdcKalHelixSegCol");
5886 if(segsc != StatusCode::SUCCESS) {
5887 log << MSG::FATAL <<
"Could not unregister RecMdcKalHelixSegCol collection" << endreq;
5892 segsc = eventSvc()->registerObject(
"/Event/Recon/RecMdcKalHelixSegCol", segcol);
5893 if( segsc.isFailure() ) {
5894 log << MSG::FATAL <<
"Could not register RecMdcKalHelixSeg" << endreq;
5897 log << MSG::INFO <<
"RecMdcKalHelixSegCol registered successfully!" <<endreq;
5900 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.);
5901 double r1(0.),r2(0.),kap1(999.),kap2(999.),tanl1(0.),tanl2(0.);
5904 SmartDataPtr<RecMdcKalTrackCol> kaltrkCol(eventSvc(),
"/Event/Recon/RecMdcKalTrackCol");
5906 log << MSG::FATAL <<
"Could not find RecMdcKalTrackCol" << endreq;
5909 log << MSG::INFO <<
"Begin to check RecMdcKalTrackCol"<<endreq;
5910 RecMdcKalTrackCol::iterator iter_trk = kaltrkCol->begin();
5911 for(
int jj=1; iter_trk != kaltrkCol->end(); iter_trk++,jj++) {
5912 log << MSG::DEBUG <<
"retrieved MDC Kalmantrack:"
5913 <<
"Track Id: " << (*iter_trk)->getTrackId()
5914 <<
" Mass of the fit: "<< (*iter_trk)->getMass(2)<< endreq
5915 <<
" Length of the track: "<< (*iter_trk)->getLength(2)
5916 <<
" Tof of the track: "<< (*iter_trk)->getTof(2) << endreq
5917 <<
" Chisq of the fit: "<< (*iter_trk)->getChisq(0,2)
5918 <<
" "<< (*iter_trk)->getChisq(1,2) << endreq
5919 <<
"Ndf of the fit: "<< (*iter_trk)->getNdf(0,1)
5920 <<
" "<< (*iter_trk)->getNdf(1,2) << endreq
5921 <<
"Kappa " << (*iter_trk)->getZHelix()[2]
5922 <<
"zhelixmu "<<(*iter_trk)->getZHelixMu()
5927 std::cout<<
"the size of gothelixsegs ..."<<gothelixsegs.size()<<std::endl;
5930 HelixSegRefVec::iterator it_gothelixseg = gothelixsegs.begin();
5931 for( ; it_gothelixseg != gothelixsegs.end(); it_gothelixseg++) {
5933 std::cout<<
"the layerId of this helixseg is ..."<<(*it_gothelixseg)->getLayerId()<<std::endl;
5934 std::cout<<
"the residual of this helixseg exclude the meas hit"<<(*it_gothelixseg)->getResExcl()<<std::endl;
5935 std::cout<<
"the residual of this helixseg include the meas hit"<<(*it_gothelixseg)->getResIncl()<<std::endl;
5936 std::cout<<
"the track id of the helixseg is ..."<<(*it_gothelixseg)->getTrackId() <<std::endl;
5937 std::cout<<
"the tof of the helixseg is ..."<<(*it_gothelixseg)->getTof()<<std::endl;
5938 std::cout<<
"the Zhit of the helixseg is ..."<<(*it_gothelixseg)->getZhit()<<std::endl;
5941 for(
int i = 0; i<43; i++) {
5942 log << MSG::DEBUG <<
"retrieved pathl["<<i<<
"]= "
5943 << (*iter_trk)->getPathl(i) <<endreq;
5947 m_trackid = (*iter_trk)->getTrackId();
5948 for(
int jj =0, iii=0; jj<5; jj++) {
5949 m_length[jj] = (*iter_trk)->getLength(jj);
5950 m_tof[jj] = (*iter_trk)->getTof(jj);
5951 m_nhits[jj] = (*iter_trk)->getNhits(jj);
5952 m_zhelix[jj] = (*iter_trk)->getZHelix()[jj];
5953 m_zhelixe[jj] = (*iter_trk)->getZHelixE()[jj];
5954 m_zhelixmu[jj] = (*iter_trk)->getZHelixMu()[jj];
5955 m_zhelixk[jj] = (*iter_trk)->getZHelixK()[jj];
5956 m_zhelixp[jj] = (*iter_trk)->getZHelixP()[jj];
5957 m_fhelix[jj] = (*iter_trk)->getFHelix()[jj];
5958 m_fhelixe[jj] = (*iter_trk)->getFHelixE()[jj];
5959 m_fhelixmu[jj] = (*iter_trk)->getFHelixMu()[jj];
5960 m_fhelixk[jj] = (*iter_trk)->getFHelixK()[jj];
5961 m_fhelixp[jj] = (*iter_trk)->getFHelixP()[jj];
5962 m_lhelix[jj] = (*iter_trk)->getLHelix()[jj];
5963 m_lhelixe[jj] = (*iter_trk)->getLHelixE()[jj];
5964 m_lhelixmu[jj] = (*iter_trk)->getLHelixMu()[jj];
5965 m_lhelixk[jj] = (*iter_trk)->getLHelixK()[jj];
5966 m_lhelixp[jj] = (*iter_trk)->getLHelixP()[jj];
5968 for(
int kk=0; kk<=jj; kk++,iii++) {
5969 m_zerror[iii] = (*iter_trk)->getZError()[jj][kk];
5970 m_zerrore[iii] = (*iter_trk)->getZErrorE()[jj][kk];
5971 m_zerrormu[iii] = (*iter_trk)->getZErrorMu()[jj][kk];
5972 m_zerrork[iii] = (*iter_trk)->getZErrorK()[jj][kk];
5973 m_zerrorp[iii] = (*iter_trk)->getZErrorP()[jj][kk];
5974 m_ferror[iii] = (*iter_trk)->getFError()[jj][kk];
5975 m_ferrore[iii] = (*iter_trk)->getFErrorE()[jj][kk];
5976 m_ferrormu[iii] = (*iter_trk)->getFErrorMu()[jj][kk];
5977 m_ferrork[iii] = (*iter_trk)->getFErrorK()[jj][kk];
5978 m_ferrorp[iii] = (*iter_trk)->getFErrorP()[jj][kk];
5979 m_lerror[iii] = (*iter_trk)->getLError()[jj][kk];
5980 m_lerrore[iii] = (*iter_trk)->getLErrorE()[jj][kk];
5981 m_lerrormu[iii] = (*iter_trk)->getLErrorMu()[jj][kk];
5982 m_lerrork[iii] = (*iter_trk)->getLErrorK()[jj][kk];
5983 m_lerrorp[iii] = (*iter_trk)->getLErrorP()[jj][kk];
5990 m_chisq[0][0] = (*iter_trk)->getChisq(0,0);
5991 m_chisq[1][0] = (*iter_trk)->getChisq(0,1);
5992 m_chisq[2][0] = (*iter_trk)->getChisq(0,2);
5993 m_chisq[3][0] = (*iter_trk)->getChisq(0,3);
5994 m_chisq[4][0] = (*iter_trk)->getChisq(0,4);
5995 m_chisq[0][1] = (*iter_trk)->getChisq(1,0);
5996 m_chisq[1][1] = (*iter_trk)->getChisq(1,1);
5997 m_chisq[2][1] = (*iter_trk)->getChisq(1,2);
5998 m_chisq[3][1] = (*iter_trk)->getChisq(1,3);
5999 m_chisq[4][1] = (*iter_trk)->getChisq(1,4);
6001 m_ndf[0][0] = (*iter_trk)->getNdf(0,0);
6002 m_ndf[1][0] = (*iter_trk)->getNdf(0,1);
6003 m_ndf[2][0] = (*iter_trk)->getNdf(0,2);
6004 m_ndf[3][0] = (*iter_trk)->getNdf(0,3);
6005 m_ndf[4][0] = (*iter_trk)->getNdf(0,4);
6006 m_ndf[0][1] = (*iter_trk)->getNdf(1,0);
6007 m_ndf[1][1] = (*iter_trk)->getNdf(1,1);
6008 m_ndf[2][1] = (*iter_trk)->getNdf(1,2);
6009 m_ndf[3][1] = (*iter_trk)->getNdf(1,3);
6010 m_ndf[4][1] = (*iter_trk)->getNdf(1,4);
6012 m_stat[0][0] = (*iter_trk)->getStat(0,0);
6013 m_stat[1][0] = (*iter_trk)->getStat(0,1);
6014 m_stat[2][0] = (*iter_trk)->getStat(0,2);
6015 m_stat[3][0] = (*iter_trk)->getStat(0,3);
6016 m_stat[4][0] = (*iter_trk)->getStat(0,4);
6017 m_stat[0][1] = (*iter_trk)->getStat(1,0);
6018 m_stat[1][1] = (*iter_trk)->getStat(1,1);
6019 m_stat[2][1] = (*iter_trk)->getStat(1,2);
6020 m_stat[3][1] = (*iter_trk)->getStat(1,3);
6021 m_stat[4][1] = (*iter_trk)->getStat(1,4);
6023 m_fptot = sqrt(1+pow(m_fhelix[4],2))/m_fhelix[2];
6024 m_fptote = sqrt(1+pow(m_fhelixe[4],2))/m_fhelixe[2];
6025 m_fptotmu = sqrt(1+pow(m_fhelixmu[4],2))/m_fhelixmu[2];
6026 m_fptotk = sqrt(1+pow(m_fhelixk[4],2))/m_fhelixk[2];
6027 m_fptotp = sqrt(1+pow(m_fhelixp[4],2))/m_fhelixp[2];
6029 m_zpt = 1/m_zhelix[2];
6030 m_zpte = 1/m_zhelixe[2];
6031 m_zptmu = 1/m_zhelixmu[2];
6032 m_zptk = 1/m_zhelixk[2];
6033 m_zptp = 1/m_zhelixp[2];
6035 m_fpt = 1/m_fhelix[2];
6036 m_fpte = 1/m_fhelixe[2];
6037 m_fptmu = 1/m_fhelixmu[2];
6038 m_fptk = 1/m_fhelixk[2];
6039 m_fptp = 1/m_fhelixp[2];
6041 m_lpt = 1/m_lhelix[2];
6042 m_lpte = 1/m_lhelixe[2];
6043 m_lptmu = 1/m_lhelixmu[2];
6044 m_lptk = 1/m_lhelixk[2];
6045 m_lptp = 1/m_lhelixp[2];
6047 m_lptot = sqrt(1+pow(m_lhelix[4],2))/m_lhelix[2];
6048 m_lptote = sqrt(1+pow(m_lhelixe[4],2))/m_lhelixe[2];
6049 m_lptotmu = sqrt(1+pow(m_lhelixmu[4],2))/m_lhelixmu[2];
6050 m_lptotk = sqrt(1+pow(m_lhelixk[4],2))/m_lhelixk[2];
6051 m_lptotp = sqrt(1+pow(m_lhelixp[4],2))/m_lhelixp[2];
6053 m_zptot = sqrt(1+pow(m_zhelix[4],2))/m_zhelix[2];
6054 m_zptote = sqrt(1+pow(m_zhelixe[4],2))/m_zhelixe[2];
6055 m_zptotmu = sqrt(1+pow(m_zhelixmu[4],2))/m_zhelixmu[2];
6056 m_zptotk = sqrt(1+pow(m_zhelixk[4],2))/m_zhelixk[2];
6057 m_zptotp = sqrt(1+pow(m_zhelixp[4],2))/m_zhelixp[2];
6059 m_zsigp = sqrt(pow((m_zptot/m_zhelix[2]),2)*m_zerror[5]+
6060 pow((m_zhelix[4]/m_zptot),2)*pow((1/m_zhelix[2]),4)*m_zerror[14]-
6061 2*m_zhelix[4]*m_zerror[12]*pow((1/m_zhelix[2]),3));
6062 m_zsigpe = sqrt(pow((m_zptote/m_zhelixe[2]),2)*m_zerrore[5]+
6063 pow((m_zhelixe[4]/m_zptote),2)*pow((1/m_zhelixe[2]),4)*m_zerrore[14]-
6064 2*m_zhelixe[4]*m_zerrore[12]*pow((1/m_zhelixe[2]),3));
6065 m_zsigpmu = sqrt(pow((m_zptotmu/m_zhelixmu[2]),2)*m_zerrormu[5]+
6066 pow((m_zhelixmu[4]/m_zptotmu),2)*pow((1/m_zhelixmu[2]),4)*m_zerrormu[14]-
6067 2*m_zhelixmu[4]*m_zerrormu[12]*pow((1/m_zhelixmu[2]),3));
6068 m_zsigpk = sqrt(pow((m_zptotk/m_zhelixk[2]),2)*m_zerrork[5]+
6069 pow((m_zhelixk[4]/m_zptotk),2)*pow((1/m_zhelixk[2]),4)*m_zerrork[14]-
6070 2*m_zhelixk[4]*m_zerrork[12]*pow((1/m_zhelixk[2]),3));
6071 m_zsigpp = sqrt(pow((m_zptotp/m_zhelixp[2]),2)*m_zerrorp[5]+
6072 pow((m_zhelixp[4]/m_zptotp),2)*pow((1/m_zhelixp[2]),4)*m_zerrorp[14]-
6073 2*m_zhelixp[4]*m_zerrorp[12]*pow((1/m_zhelixp[2]),3));
6076 StatusCode sc1 = m_nt1->write();
6077 if( sc1.isFailure() ) cout<<
"Ntuple1 filling failed!"<<endl;
6082 phi1 = (*iter_trk)->getFFi0();
6083 r1 = (*iter_trk)->getFDr();
6084 z1 = (*iter_trk)->getFDz();
6085 kap1 = (*iter_trk)->getFCpa();
6086 tanl1 = (*iter_trk)->getFTanl();
6089 p1 = sqrt(1+tanl1*tanl1)/kap1;
6090 the1 =
M_PI/2-atan(tanl1);
6091 }
else if(jj == 2) {
6092 phi2 = (*iter_trk)->getFFi0();
6093 r2 = (*iter_trk)->getFDr();
6094 z2 = (*iter_trk)->getFDz();
6095 kap2 = (*iter_trk)->getFCpa();
6096 tanl2 = (*iter_trk)->getFTanl();
6099 p2 = sqrt(1+tanl2*tanl2)/kap1;
6100 the2 =
M_PI/2-atan(tanl2);
6108 m_delthe = the1 + the2;
6111 StatusCode sc2 = m_nt2->write();
6112 if( sc2.isFailure() ) cout<<
"Ntuple2 filling failed!"<<endl;
6120 cout <<
"Kalfitting finished " << std::endl;
7381 std::map<double,double> esti1_r;
7382 std::map<double,double> esti2_r;
7383 std::map<double,double> esti1_phi;
7384 std::map<double,double> esti2_phi;
7385 std::map<double,double> esti1_z;
7386 std::map<double,double> esti2_z;
7387 std::map<double,double> diff1_phi;
7388 std::map<double,double> diff1_z;
7389 std::map<double,double> diff2_phi;
7390 std::map<double,double> diff2_z;
7391 std::map<double,HepPoint3D> posnew;
7392 std::map<double,HepSymMatrix> eanew;
7393 std::map<double,HepVector> anew;
7394 std::map<double,double> meas_r;
7395 std::map<double,double> meas_phi;
7396 std::map<double,double> meas_z;
7397 std::map<double,double> meas_phierr;
7398 std::map<double,double> esti_phierr;
7399 std::map<double,double> meas_zerr;
7400 std::map<double,double> esti_zerr;
7402 if(
debug_==4){ cout <<
"wall size" << _BesKalmanFitWalls.size() << endl;}
7403 for(
int iLayer=2; iLayer>=0; iLayer--)
7421 meas_phierr.clear();
7423 esti_phierr.clear();
7426 int nHits=myGemHitCol[iLayer].size();
7429 if (nHits <= 0)
continue;
7432 vector<KalFitGemHit>::iterator
iter;
7436 innerwall(track,hypo,way,64+iLayer*(-32),88+iLayer*(-32));
7439 for(
iter=myGemHitCol[iLayer].begin();
iter!=myGemHitCol[iLayer].end();
iter++ )
7442 track.
set(track_start.
pivot(),track_start.
a(),track_start.
Ea());
7444 double rGem=(*iter).getR();
7445 double phiGem=(*iter).getPhi();
7446 double zGem=(*iter).getZ();
7447 HepVector v_measu=(*iter).getVecPos();
7460 double x0=posEstiAtGem.x();
7461 double y0=posEstiAtGem.y();
7469 if(
muls_) track.
ms(pathInGem, _BesKalmanFitMaterials[26+32*iLayer], way);
7470 if(
loss_) track.
eloss(pathInGem, _BesKalmanFitMaterials[26+32*iLayer], way);
7473 double rGemEsti=posEstiAtGem.perp();
7474 double phiGemEsti=posEstiAtGem.phi();
7475 double zGemEsti=posEstiAtGem.z();
7476 HepVector v_estim(2,0);
7477 v_estim(1)=phiGemEsti;
7478 v_estim(2)=zGemEsti;
7485 const HepSymMatrix& Ea = track.
Ea();
7486 HepVector v_a = track.
a();
7489 double kappa=v_a(3);
7494 HepMatrix
H(2, 5, 0);
7497 H(1,1)=-y0*
cos(phi0)/(y0*y0+x0*x0)+x0*
sin(phi0)/(x0*x0+y0*y0);
7500 H(1,2)=-y0/(y0*y0+x0*x0)*((-1)*drho*
sin(phi0)+
m_alpha/kappa*(
sin(phi0+dPhi)-
sin(phi0)))+x0/(x0*x0+y0*y0)*(drho*
cos(phi0)+
m_alpha/kappa*(
cos(phi0)-
cos(phi0+dPhi)));
7501 H(1,3)=-y0/(y0*y0+x0*x0)*(-1)*
m_alpha/(kappa*kappa)*(
cos(phi0)-
cos(phi0+dPhi))+x0/(x0*x0+y0*y0)*(-1)*
m_alpha/(kappa*kappa)*(
sin(phi0)-
sin(phi0+dPhi));
7502 H(2,3)=
m_alpha/(kappa*kappa)*tl*dPhi;
7506 HepMatrix H_T=
H.T();
7511 HepSymMatrix V=(*iter).getErrMat();
7514 HepMatrix HEaH_T=HEa*H_T;
7516 HepMatrix Vinv=(V+HEaH_T).inverse(ierr);
7517 if(ierr!=0) cout<<
"error in HEaH_T.inverse()!"<<endl;
7519 HepMatrix K=Ea*H_T*Vinv;
7521 HepSymMatrix EaNew(5,0);
7522 EaNew.assign(Ea-K*
H*Ea);
7524 HepVector v_diff=v_measu-v_estim;
7526 double delPhi=v_diff(1);
7527 if(fabs(v_diff(1))>6.28) v_measu(1)=-1*v_measu(1);
7534 HepVector aNew=v_a+K*v_diff;
7540 HepPoint3D posEstiAtGem_new = track_test.
x(dPhiToGem_New);
7541 double phiGemEsti_new=posEstiAtGem_new.phi();
7542 double zGemEsti_new=posEstiAtGem_new.z();
7543 double x0_new = posEstiAtGem_new.x();
7544 double y0_new = posEstiAtGem_new.y();
7545 HepVector v_estim_new(2,0);
7546 v_estim_new(1)=phiGemEsti_new;
7547 v_estim_new(2)=zGemEsti_new;
7551 v_diff=v_measu-v_estim_new;
7559 track_test.
pivot_numf(posEstiAtGem_new, pathInGem);
7560 HepVector v_a_new = track_test.
a();
7562 HepSymMatrix Ea_new = track_test.
Ea();
7564 double phi0_new = v_a_new(2);
7565 double kappa_new=v_a_new(3);
7566 HepMatrix H_new(2, 5, 0);
7575 HepMatrix H_new_T=H_new.T();
7582 HepSymMatrix R(2,0);
7584 R.assign(V-
H*EaNew*H_T);
7592 HepVector v_dChi2=v_diff.T()*R.inverse(ierr)*v_diff;
7593 if(ierr!=0) cout<<
"error in R.inverse()!"<<endl;
7596 double dChi2=v_dChi2(1);
7600 cout<<
"pivot: "<<posEstiAtGem<<endl;
7601 cout<<
"new pivot: "<<posEstiAtGem<<endl;
7602 cout<<
"a: "<<v_a<<endl;
7603 cout<<
"new a: "<<aNew<<endl;
7604 cout<<
"Ea: "<<Ea<<endl;
7605 cout<<
"new Ea: "<<EaNew<<endl;
7606 cout<<
"dchi2= "<<dChi2<<endl;
7641 esti1_r[dChi2]=rGemEsti;
7642 esti1_phi[dChi2]=phiGemEsti;
7643 esti1_z[dChi2]=zGemEsti;
7644 esti2_r[dChi2]=rGem;
7645 esti2_phi[dChi2]=phiGemEsti_new;
7646 esti2_z[dChi2]=zGemEsti_new;
7647 diff1_phi[dChi2]=phiGem-phiGemEsti;
7648 diff1_z[dChi2]=zGem-zGemEsti;
7649 diff2_phi[dChi2]=phiGem-phiGemEsti_new;
7650 diff2_z[dChi2]=zGem-zGemEsti_new;
7651 posnew[dChi2]=posEstiAtGem;
7655 meas_phi[dChi2]=phiGem;
7658 meas_phierr[dChi2]=sqrt(V[0][0]);
7660 meas_zerr[dChi2]=sqrt(V[1][1]);
7662 esti_phierr[dChi2]=sqrt((
H*(track_start.
Ea())*H_T)[0][0]);
7664 esti_zerr[dChi2]=sqrt((
H*(track_start.
Ea())*H_T)[1][1]);
7673 track.
set((posnew.begin())->second,(anew.begin())->second,(eanew.begin()->second));
7677 HepPoint3D new_posEstiAtGem = track.
x(new_dPhiToGem);
7678 double new_pathInGem;
7681 track.
pivot_numf(new_posEstiAtGem, new_pathInGem);
7686 if(
muls_) track.
ms(new_pathInGem, _BesKalmanFitMaterials[26+32*iLayer], way);
7687 if(
loss_) track.
eloss(new_pathInGem, _BesKalmanFitMaterials[26+32*iLayer], way);
7691 innerwall(track,hypo,way,90+(-32)*iLayer,95+(-32)*iLayer);}
7693 innerwall(track,hypo,way,90+(-32)*iLayer,94+(-32)*iLayer);}
7701 m_esti1_r[ii]=(esti1_r.begin())->second;
7703 m_esti1_phi[ii]=(esti1_phi.begin())->second;
7704 m_esti1_z[ii]=(esti1_z.begin())->second;
7705 m_esti2_r[ii]=(esti2_r.begin())->second;
7706 m_esti2_phi[ii]=(esti2_phi.begin())->second;
7707 m_esti2_z[ii]=(esti2_z.begin())->second;
7708 m_diff1_phi[ii]=(diff1_phi.begin())->second;
7709 m_diff1_z[ii]=(diff1_z.begin())->second;
7710 m_diff2_phi[ii]=(diff2_phi.begin())->second;
7711 m_diff2_z[ii]=(diff2_z.begin())->second;
7712 m_Gchi2[ii]=(esti1_r.begin())->first;
7713 m_meas_r[ii]=(meas_r.begin())->second;
7714 m_meas_phi[ii]=(meas_phi.begin())->second;
7716 m_meas_z[ii]=(meas_z.begin())->second;
7717 m_meas_phierr[ii]=(meas_phierr.begin())->second;
7718 m_meas_zerr[ii]=(meas_zerr.begin())->second;
7719 m_esti_phierr[ii]=(esti_phierr.begin())->second;
7720 m_esti_zerr[ii]=(esti_zerr.begin())->second;
7721 m_GemLayer[ii]=iLayer;
7733 if(m_meas_phi[0]>0 && m_meas_phi[0]<1.5707963) m_qua=0;
7734 if(m_meas_phi[0]>1.5707963 && m_meas_phi[0]<3.1415926) m_qua=1;
7735 if(m_meas_phi[0]>-3.1415926&& m_meas_phi[0]<-1.5707963) m_qua=2;
7736 if(m_meas_phi[0]>-1.5707963&& m_meas_phi[0]<0) m_qua=3;
7759 meas_phierr.clear();
7761 esti_phierr.clear();
7767 SmartDataPtr<Event::EventHeader> evtHeader(eventSvc(),
"/Event/EventHeader");
7769 m_evt3 = evtHeader->eventNumber();
7781 SmartDataPtr<Event::EventHeader> eventHeader(eventSvc(),
"/Event/EventHeader");
7783 cout <<
"Could not find Event Header" << endl;
7785 int eventNo = eventHeader->eventNumber();
7786 int runNo = eventHeader->runNumber();
7787 SmartDataPtr<RecMdcTrackCol> recMdcTrack(eventSvc(),
"/Event/Recon/RecMdcTrackCol");
7789 cout <<
"Could not find RecMdcTrackCol" << endl;
7791 RecMdcTrackCol::iterator itrk = recMdcTrack->begin();
7793 for(; itrk != recMdcTrack->end(); itrk++){
7794 if((*itrk)->trackId() == track.
trasan_id() && (*itrk)->getNcluster() > 0) {
7795 clusterRefVec = (*itrk)->getVecClusters();
7803 if(
debug_==4){ cout<<
"wall size"<<_BesKalmanFitWalls.size()<<endl;}
7806 ClusterRefVec::iterator itCluster;
7807 ClusterRefVec::iterator itStartFlag;
7808 ClusterRefVec::iterator itStopFlag;
7812 double R_o_cgem=_BesKalmanFitWalls[0].radius();
7814 itStartFlag = clusterRefVec.begin();
7815 itStopFlag = clusterRefVec.end();
7820 itStartFlag = clusterRefVec.end()-1;
7821 itStopFlag = clusterRefVec.begin()-1;
7831 for(itCluster = itStartFlag; itCluster !=itStopFlag; itCluster = itCluster + step){
7832 int layer = (*itCluster)->getlayerid();
7835 innerwall(track,hypo,way,lastR,dmR,index);
7845 m_dchi2_hit_e[m_hit_no] = dchi2;
7846 m_residest_hit_e[m_hit_no] = 0;
7847 m_residnew_hit_e[m_hit_no] = 0;
7848 m_layer_hit_e[m_hit_no] = layer;
7850 m_anal_dr_hit_e[m_hit_no] = worktemp.
a()[0];
7851 m_anal_phi0_hit_e[m_hit_no] = worktemp.
a()[1];
7852 m_anal_kappa_hit_e[m_hit_no] = worktemp.
a()[2];
7853 m_anal_dz_hit_e[m_hit_no] = worktemp.
a()[3];
7854 m_anal_tanl_hit_e[m_hit_no] = worktemp.
a()[4];
7855 m_anal_ea_dr_hit_e[m_hit_no] = worktemp.
Ea()[0][0];
7856 m_anal_ea_phi0_hit_e[m_hit_no] = worktemp.
Ea()[1][1];
7857 m_anal_ea_kappa_hit_e[m_hit_no] = worktemp.
Ea()[2][2];
7858 m_anal_ea_dz_hit_e[m_hit_no] = worktemp.
Ea()[3][3];
7859 m_anal_ea_tanl_hit_e[m_hit_no] = worktemp.
Ea()[4][4];
7864 m_dchi2_hit_mu[m_hit_no] = dchi2;
7865 m_residest_hit_mu[m_hit_no] = 0;
7866 m_residnew_hit_mu[m_hit_no] = 0;
7867 m_layer_hit_mu[m_hit_no] = layer;
7869 m_anal_dr_hit_mu[m_hit_no] = worktemp.
a()[0];
7870 m_anal_phi0_hit_mu[m_hit_no] = worktemp.
a()[1];
7871 m_anal_kappa_hit_mu[m_hit_no] = worktemp.
a()[2];
7872 m_anal_dz_hit_mu[m_hit_no] = worktemp.
a()[3];
7873 m_anal_tanl_hit_mu[m_hit_no] = worktemp.
a()[4];
7874 m_anal_ea_dr_hit_mu[m_hit_no] = worktemp.
Ea()[0][0];
7875 m_anal_ea_phi0_hit_mu[m_hit_no] = worktemp.
Ea()[1][1];
7876 m_anal_ea_kappa_hit_mu[m_hit_no] = worktemp.
Ea()[2][2];
7877 m_anal_ea_dz_hit_mu[m_hit_no] = worktemp.
Ea()[3][3];
7878 m_anal_ea_tanl_hit_mu[m_hit_no] = worktemp.
Ea()[4][4];
7884 m_dchi2_hit_pi[m_hit_no] = dchi2;
7885 m_residest_hit_pi[m_hit_no] = 0;
7886 m_residnew_hit_pi[m_hit_no] = 0;
7887 m_layer_hit_pi[m_hit_no] = layer;
7889 m_anal_dr_hit_pi[m_hit_no] = worktemp.
a()[0];
7890 m_anal_phi0_hit_pi[m_hit_no] = worktemp.
a()[1];
7891 m_anal_kappa_hit_pi[m_hit_no] = worktemp.
a()[2];
7892 m_anal_dz_hit_pi[m_hit_no] = worktemp.
a()[3];
7893 m_anal_tanl_hit_pi[m_hit_no] = worktemp.
a()[4];
7894 m_anal_ea_dr_hit_pi[m_hit_no] = worktemp.
Ea()[0][0];
7895 m_anal_ea_phi0_hit_pi[m_hit_no] = worktemp.
Ea()[1][1];
7896 m_anal_ea_kappa_hit_pi[m_hit_no] = worktemp.
Ea()[2][2];
7897 m_anal_ea_dz_hit_pi[m_hit_no] = worktemp.
Ea()[3][3];
7898 m_anal_ea_tanl_hit_pi[m_hit_no] = worktemp.
Ea()[4][4];
7903 m_dchi2_hit_k[m_hit_no] = dchi2;
7904 m_residest_hit_k[m_hit_no] = 0;
7905 m_residnew_hit_k[m_hit_no] = 0;
7906 m_layer_hit_k[m_hit_no] = layer;
7908 m_anal_dr_hit_k[m_hit_no] = worktemp.
a()[0];
7909 m_anal_phi0_hit_k[m_hit_no] = worktemp.
a()[1];
7910 m_anal_kappa_hit_k[m_hit_no] = worktemp.
a()[2];
7911 m_anal_dz_hit_k[m_hit_no] = worktemp.
a()[3];
7912 m_anal_tanl_hit_k[m_hit_no] = worktemp.
a()[4];
7913 m_anal_ea_dr_hit_k[m_hit_no] = worktemp.
Ea()[0][0];
7914 m_anal_ea_phi0_hit_k[m_hit_no] = worktemp.
Ea()[1][1];
7915 m_anal_ea_kappa_hit_k[m_hit_no] = worktemp.
Ea()[2][2];
7916 m_anal_ea_dz_hit_k[m_hit_no] = worktemp.
Ea()[3][3];
7917 m_anal_ea_tanl_hit_k[m_hit_no] = worktemp.
Ea()[4][4];
7922 m_dchi2_hit_p[m_hit_no] = dchi2;
7923 m_residest_hit_p[m_hit_no] = 0;
7924 m_residnew_hit_p[m_hit_no] = 0;
7925 m_layer_hit_p[m_hit_no] = layer;
7927 m_anal_dr_hit_p[m_hit_no] = worktemp.
a()[0];
7928 m_anal_phi0_hit_p[m_hit_no] = worktemp.
a()[1];
7929 m_anal_kappa_hit_p[m_hit_no] = worktemp.
a()[2];
7930 m_anal_dz_hit_p[m_hit_no] = worktemp.
a()[3];
7931 m_anal_tanl_hit_p[m_hit_no] = worktemp.
a()[4];
7932 m_anal_ea_dr_hit_p[m_hit_no] = worktemp.
Ea()[0][0];
7933 m_anal_ea_phi0_hit_p[m_hit_no] = worktemp.
Ea()[1][1];
7934 m_anal_ea_kappa_hit_p[m_hit_no] = worktemp.
Ea()[2][2];
7935 m_anal_ea_dz_hit_p[m_hit_no] = worktemp.
Ea()[3][3];
7936 m_anal_ea_tanl_hit_p[m_hit_no] = worktemp.
Ea()[4][4];
7944 innerwall(track,hypo,way,lastR,R_o_cgem,index);
7948 innerwall(track,hypo,way,lastR,0,index);