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(
"mchelix",5,m_mchelix);
581 status7 = m_nt7->addItem(
"mcptot",m_mcptot);
582 status7 = m_nt7->addItem(
"mcpid",m_mcpid);
583 status7 = m_nt7 ->addItem(
"dropedHits_kal_e",m_dropedHits_kal_e);
584 status7 = m_nt7 ->addItem(
"kappa2_kal_e",m_kappa2_kal_e);
585 status7 = m_nt7 ->addItem(
"trackNhits_kal_e",m_trackNhits_kal_e);
586 status7 = m_nt7 ->addItem(
"trackNster_kal_e",m_trackNster_kal_e);
587 status7 = m_nt7 ->addItem(
"trackNaxis_kal_e",m_trackNaxis_kal_e);
588 status7 = m_nt7 ->addItem(
"chi2_kal_e",m_chi2_kal_e);
589 status7 = m_nt7 ->addItem(
"Ea00_kal_e",m_Ea00_kal_e);
590 status7 = m_nt7 ->addItem(
"Ea11_kal_e",m_Ea11_kal_e);
591 status7 = m_nt7 ->addItem(
"Ea22_kal_e",m_Ea22_kal_e);
592 status7 = m_nt7 ->addItem(
"Ea33_kal_e",m_Ea33_kal_e);
593 status7 = m_nt7 ->addItem(
"Ea44_kal_e",m_Ea44_kal_e);
594 status7 = m_nt7 ->addItem(
"dropedHits_kal_mu",m_dropedHits_kal_mu);
595 status7 = m_nt7 ->addItem(
"kappa2_kal_mu",m_kappa2_kal_mu);
596 status7 = m_nt7 ->addItem(
"trackNhits_kal_mu",m_trackNhits_kal_mu);
597 status7 = m_nt7 ->addItem(
"trackNster_kal_mu",m_trackNster_kal_mu);
598 status7 = m_nt7 ->addItem(
"trackNaxis_kal_mu",m_trackNaxis_kal_mu);
599 status7 = m_nt7 ->addItem(
"chi2_kal_mu",m_chi2_kal_mu);
600 status7 = m_nt7 ->addItem(
"Ea00_kal_mu",m_Ea00_kal_mu);
601 status7 = m_nt7 ->addItem(
"Ea11_kal_mu",m_Ea11_kal_mu);
602 status7 = m_nt7 ->addItem(
"Ea22_kal_mu",m_Ea22_kal_mu);
603 status7 = m_nt7 ->addItem(
"Ea33_kal_mu",m_Ea33_kal_mu);
604 status7 = m_nt7 ->addItem(
"Ea44_kal_mu",m_Ea44_kal_mu);
605 status7 = m_nt7 ->addItem(
"iqual_front_kal_mu",m_iqual_front_kal_mu);
606 status7 = m_nt7 ->addItem(
"dropedHits_kal_pi",m_dropedHits_kal_pi);
607 status7 = m_nt7 ->addItem(
"kappa2_kal_pi",m_kappa2_kal_pi);
608 status7 = m_nt7 ->addItem(
"trackNhits_kal_pi",m_trackNhits_kal_pi);
609 status7 = m_nt7 ->addItem(
"trackNster_kal_pi",m_trackNster_kal_pi);
610 status7 = m_nt7 ->addItem(
"trackNaxis_kal_pi",m_trackNaxis_kal_pi);
611 status7 = m_nt7 ->addItem(
"chi2_kal_pi",m_chi2_kal_pi);
612 status7 = m_nt7 ->addItem(
"Ea00_kal_pi",m_Ea00_kal_pi);
613 status7 = m_nt7 ->addItem(
"Ea11_kal_pi",m_Ea11_kal_pi);
614 status7 = m_nt7 ->addItem(
"Ea22_kal_pi",m_Ea22_kal_pi);
615 status7 = m_nt7 ->addItem(
"Ea33_kal_pi",m_Ea33_kal_pi);
616 status7 = m_nt7 ->addItem(
"Ea44_kal_pi",m_Ea44_kal_pi);
617 status7 = m_nt7 ->addItem(
"dropedHits_kal_k",m_dropedHits_kal_k);
618 status7 = m_nt7 ->addItem(
"kappa2_kal_k",m_kappa2_kal_k);
619 status7 = m_nt7 ->addItem(
"trackNhits_kal_k",m_trackNhits_kal_k);
620 status7 = m_nt7 ->addItem(
"trackNster_kal_k",m_trackNster_kal_k);
621 status7 = m_nt7 ->addItem(
"trackNaxis_kal_k",m_trackNaxis_kal_k);
622 status7 = m_nt7 ->addItem(
"chi2_kal_k",m_chi2_kal_k);
623 status7 = m_nt7 ->addItem(
"Ea00_kal_k",m_Ea00_kal_k);
624 status7 = m_nt7 ->addItem(
"Ea11_kal_k",m_Ea11_kal_k);
625 status7 = m_nt7 ->addItem(
"Ea22_kal_k",m_Ea22_kal_k);
626 status7 = m_nt7 ->addItem(
"Ea33_kal_k",m_Ea33_kal_k);
627 status7 = m_nt7 ->addItem(
"Ea44_kal_k",m_Ea44_kal_k);
628 status7 = m_nt7 ->addItem(
"dropedHits_kal_p",m_dropedHits_kal_p);
629 status7 = m_nt7 ->addItem(
"kappa2_kal_p",m_kappa2_kal_p);
630 status7 = m_nt7 ->addItem(
"trackNhits_kal_p",m_trackNhits_kal_p);
631 status7 = m_nt7 ->addItem(
"trackNster_kal_p",m_trackNster_kal_p);
632 status7 = m_nt7 ->addItem(
"trackNaxis_kal_p",m_trackNaxis_kal_p);
633 status7 = m_nt7 ->addItem(
"chi2_kal_p",m_chi2_kal_p);
634 status7 = m_nt7 ->addItem(
"Ea00_kal_p",m_Ea00_kal_p);
635 status7 = m_nt7 ->addItem(
"Ea11_kal_p",m_Ea11_kal_p);
636 status7 = m_nt7 ->addItem(
"Ea22_kal_p",m_Ea22_kal_p);
637 status7 = m_nt7 ->addItem(
"Ea33_kal_p",m_Ea33_kal_p);
638 status7 = m_nt7 ->addItem(
"Ea44_kal_p",m_Ea44_kal_p);
640 status7 = m_nt7 ->addItem(
"hit_number",m_hit_no,0, 1000);
641 status7 = m_nt7 ->addItem(
"nCluster",m_nCluster,0, 1000);
642 status7 = m_nt7->addIndexedItem(
"dchi2_hit_e",m_hit_no,m_dchi2_hit_e);
643 status7 = m_nt7->addIndexedItem(
"residual_estim_hit_e",m_hit_no,m_residest_hit_e);
644 status7 = m_nt7->addIndexedItem(
"residual_hit_e",m_hit_no,m_residnew_hit_e);
645 status7 = m_nt7->addIndexedItem(
"layer_hit_e",m_hit_no,m_layer_hit_e);
646 status7 = m_nt7->addIndexedItem(
"kaldr_hit_e",m_hit_no,m_anal_dr_hit_e);
647 status7 = m_nt7->addIndexedItem(
"kalphi0_hit_e",m_hit_no,m_anal_phi0_hit_e);
648 status7 = m_nt7->addIndexedItem(
"kalkappa_hit_e",m_hit_no,m_anal_kappa_hit_e);
649 status7 = m_nt7->addIndexedItem(
"kaldz_hit_e",m_hit_no,m_anal_dz_hit_e);
650 status7 = m_nt7->addIndexedItem(
"kaltanl_hit_e",m_hit_no,m_anal_tanl_hit_e);
651 status7 = m_nt7->addIndexedItem(
"dr_ea_hit_e",m_hit_no,m_anal_ea_dr_hit_e);
652 status7 = m_nt7->addIndexedItem(
"phi0_ea_hit_e",m_hit_no,m_anal_ea_phi0_hit_e);
653 status7 = m_nt7->addIndexedItem(
"kappa_ea_hit_e",m_hit_no,m_anal_ea_kappa_hit_e);
654 status7 = m_nt7->addIndexedItem(
"dz_ea_hit_e",m_hit_no,m_anal_ea_dz_hit_e);
655 status7 = m_nt7->addIndexedItem(
"tanl_ea_hit_e",m_hit_no,m_anal_ea_tanl_hit_e);
656 status7 = m_nt7->addIndexedItem(
"dchi2_hit_mu",m_hit_no,m_dchi2_hit_mu);
657 status7 = m_nt7->addIndexedItem(
"residual_estim_hit_mu",m_hit_no,m_residest_hit_mu);
658 status7 = m_nt7->addIndexedItem(
"residual_hit_mu",m_hit_no,m_residnew_hit_mu);
659 status7 = m_nt7->addIndexedItem(
"layer_hit_mu",m_hit_no,m_layer_hit_mu);
660 status7 = m_nt7->addIndexedItem(
"kaldr_hit_mu",m_hit_no,m_anal_dr_hit_mu);
661 status7 = m_nt7->addIndexedItem(
"kalphi0_hit_mu",m_hit_no,m_anal_phi0_hit_mu);
662 status7 = m_nt7->addIndexedItem(
"kalkappa_hit_mu",m_hit_no,m_anal_kappa_hit_mu);
663 status7 = m_nt7->addIndexedItem(
"kaldz_hit_mu",m_hit_no,m_anal_dz_hit_mu);
664 status7 = m_nt7->addIndexedItem(
"kaltanl_hit_mu",m_hit_no,m_anal_tanl_hit_mu);
665 status7 = m_nt7->addIndexedItem(
"dr_ea_hit_mu",m_hit_no,m_anal_ea_dr_hit_mu);
666 status7 = m_nt7->addIndexedItem(
"phi0_ea_hit_mu",m_hit_no,m_anal_ea_phi0_hit_mu);
667 status7 = m_nt7->addIndexedItem(
"kappa_ea_hit_mu",m_hit_no,m_anal_ea_kappa_hit_mu);
668 status7 = m_nt7->addIndexedItem(
"dz_ea_hit_mu",m_hit_no,m_anal_ea_dz_hit_mu);
669 status7 = m_nt7->addIndexedItem(
"tanl_ea_hit_mu",m_hit_no,m_anal_ea_tanl_hit_mu);
670 status7 = m_nt7->addIndexedItem(
"dchi2_hit_pi",m_hit_no,m_dchi2_hit_pi);
671 status7 = m_nt7->addIndexedItem(
"residual_estim_hit_pi",m_hit_no,m_residest_hit_pi);
672 status7 = m_nt7->addIndexedItem(
"residual_hit_pi",m_hit_no,m_residnew_hit_pi);
673 status7 = m_nt7->addIndexedItem(
"layer_hit_pi",m_hit_no,m_layer_hit_pi);
674 status7 = m_nt7->addIndexedItem(
"kaldr_hit_pi",m_hit_no,m_anal_dr_hit_pi);
675 status7 = m_nt7->addIndexedItem(
"kalphi0_hit_pi",m_hit_no,m_anal_phi0_hit_pi);
676 status7 = m_nt7->addIndexedItem(
"kalkappa_hit_pi",m_hit_no,m_anal_kappa_hit_pi);
677 status7 = m_nt7->addIndexedItem(
"kaldz_hit_pi",m_hit_no,m_anal_dz_hit_pi);
678 status7 = m_nt7->addIndexedItem(
"kaltanl_hit_pi",m_hit_no,m_anal_tanl_hit_pi);
679 status7 = m_nt7->addIndexedItem(
"dr_ea_hit_pi",m_hit_no,m_anal_ea_dr_hit_pi);
680 status7 = m_nt7->addIndexedItem(
"phi0_ea_hit_pi",m_hit_no,m_anal_ea_phi0_hit_pi);
681 status7 = m_nt7->addIndexedItem(
"kappa_ea_hit_pi",m_hit_no,m_anal_ea_kappa_hit_pi);
682 status7 = m_nt7->addIndexedItem(
"dz_ea_hit_pi",m_hit_no,m_anal_ea_dz_hit_pi);
683 status7 = m_nt7->addIndexedItem(
"tanl_ea_hit_pi",m_hit_no,m_anal_ea_tanl_hit_pi);
684 status7 = m_nt7->addIndexedItem(
"dchi2_hit_k",m_hit_no,m_dchi2_hit_k);
685 status7 = m_nt7->addIndexedItem(
"residual_estim_hit_k",m_hit_no,m_residest_hit_k);
686 status7 = m_nt7->addIndexedItem(
"residual_hit_k",m_hit_no,m_residnew_hit_k);
687 status7 = m_nt7->addIndexedItem(
"layer_hit_k",m_hit_no,m_layer_hit_k);
688 status7 = m_nt7->addIndexedItem(
"kaldr_hit_k",m_hit_no,m_anal_dr_hit_k);
689 status7 = m_nt7->addIndexedItem(
"kalphi0_hit_k",m_hit_no,m_anal_phi0_hit_k);
690 status7 = m_nt7->addIndexedItem(
"kalkappa_hit_k",m_hit_no,m_anal_kappa_hit_k);
691 status7 = m_nt7->addIndexedItem(
"kaldz_hit_k",m_hit_no,m_anal_dz_hit_k);
692 status7 = m_nt7->addIndexedItem(
"kaltanl_hit_k",m_hit_no,m_anal_tanl_hit_k);
693 status7 = m_nt7->addIndexedItem(
"dr_ea_hit_k",m_hit_no,m_anal_ea_dr_hit_k);
694 status7 = m_nt7->addIndexedItem(
"phi0_ea_hit_k",m_hit_no,m_anal_ea_phi0_hit_k);
695 status7 = m_nt7->addIndexedItem(
"kappa_ea_hit_k",m_hit_no,m_anal_ea_kappa_hit_k);
696 status7 = m_nt7->addIndexedItem(
"dz_ea_hit_k",m_hit_no,m_anal_ea_dz_hit_k);
697 status7 = m_nt7->addIndexedItem(
"tanl_ea_hit_k",m_hit_no,m_anal_ea_tanl_hit_k);
698 status7 = m_nt7->addIndexedItem(
"dchi2_hit_p",m_hit_no,m_dchi2_hit_p);
699 status7 = m_nt7->addIndexedItem(
"residual_estim_hit_p",m_hit_no,m_residest_hit_p);
700 status7 = m_nt7->addIndexedItem(
"residual_hit_p",m_hit_no,m_residnew_hit_p);
701 status7 = m_nt7->addIndexedItem(
"layer_hit_p",m_hit_no,m_layer_hit_p);
702 status7 = m_nt7->addIndexedItem(
"kaldr_hit_p",m_hit_no,m_anal_dr_hit_p);
703 status7 = m_nt7->addIndexedItem(
"kalphi0_hit_p",m_hit_no,m_anal_phi0_hit_p);
704 status7 = m_nt7->addIndexedItem(
"kalkappa_hit_p",m_hit_no,m_anal_kappa_hit_p);
705 status7 = m_nt7->addIndexedItem(
"kaldz_hit_p",m_hit_no,m_anal_dz_hit_p);
706 status7 = m_nt7->addIndexedItem(
"kaltanl_hit_p",m_hit_no,m_anal_tanl_hit_p);
707 status7 = m_nt7->addIndexedItem(
"dr_ea_hit_p",m_hit_no,m_anal_ea_dr_hit_p);
708 status7 = m_nt7->addIndexedItem(
"phi0_ea_hit_p",m_hit_no,m_anal_ea_phi0_hit_p);
709 status7 = m_nt7->addIndexedItem(
"kappa_ea_hit_p",m_hit_no,m_anal_ea_kappa_hit_p);
710 status7 = m_nt7->addIndexedItem(
"dz_ea_hit_p",m_hit_no,m_anal_ea_dz_hit_p);
711 status7 = m_nt7->addIndexedItem(
"tanl_ea_hit_p",m_hit_no,m_anal_ea_tanl_hit_p);
713 if( status7.isFailure() ) cout<<
"Ntuple7 add item failed!"<<endl;
720 NTuplePtr nt6(
ntupleSvc(),
"FILE104/n106");
722 if ( nt6 ) m_nt6 = nt6;
724 m_nt6=
ntupleSvc()->book(
"FILE104/n106",CLID_ColumnWiseTuple,
"kal seg");
726 status6 = m_nt6->addItem(
"docaInc",m_docaInc);
727 status6 = m_nt6->addItem(
"docaExc",m_docaExc);
728 status6 = m_nt6->addItem(
"tdr",m_tdrift);
729 status6 = m_nt6->addItem(
"layerid", m_layerid);
730 status6 = m_nt6->addItem(
"event", m_eventNo);
731 status6 = m_nt6->addItem(
"residualInc",m_residualInc);
732 status6 = m_nt6->addItem(
"residualExc",m_residualExc);
733 status6 = m_nt6->addItem(
"lr",m_lr);
734 status6 = m_nt6->addItem(
"dd",m_dd);
735 status6 = m_nt6->addItem(
"y",m_yposition);
737 if( status6.isFailure() ) cout<<
"Ntuple6 add item failed!"<<endl;
742 NTuplePtr nt10(
ntupleSvc(),
"FILE104/n110");
744 if ( nt10 ) m_nt10 = nt10;
746 m_nt10=
ntupleSvc()->book(
"FILE104/n110",CLID_ColumnWiseTuple,
"test");
748 status10 = m_nt10->addItem(
"evt",m_evt3);
749 status10 = m_nt10->addItem(
"qua",m_qua);
750 status10 = m_nt10->addItem(
"nhit",m_nGemHits,0,4000000);
751 status10 = m_nt10->addIndexedItem(
"meas_r",m_nGemHits,m_meas_r);
752 status10 = m_nt10->addIndexedItem(
"meas_phi",m_nGemHits,m_meas_phi);
753 status10 = m_nt10->addIndexedItem(
"meas_z",m_nGemHits,m_meas_z);
754 status10 = m_nt10->addIndexedItem(
"esti1_r",m_nGemHits,m_esti1_r);
755 status10 = m_nt10->addIndexedItem(
"esti1_phi",m_nGemHits,m_esti1_phi);
756 status10 = m_nt10->addIndexedItem(
"esti1_z",m_nGemHits,m_esti1_z);
757 status10 = m_nt10->addIndexedItem(
"esti2_r",m_nGemHits,m_esti2_r);
758 status10 = m_nt10->addIndexedItem(
"esti2_phi",m_nGemHits,m_esti2_phi);
759 status10 = m_nt10->addIndexedItem(
"esti2_z",m_nGemHits,m_esti2_z);
760 status10 = m_nt10->addIndexedItem(
"diff1_phi",m_nGemHits,m_diff1_phi);
761 status10 = m_nt10->addIndexedItem(
"diff1_z",m_nGemHits,m_diff1_z);
762 status10 = m_nt10->addIndexedItem(
"diff2_phi",m_nGemHits,m_diff2_phi);
763 status10 = m_nt10->addIndexedItem(
"diff2_z",m_nGemHits,m_diff2_z);
764 status10 = m_nt10->addIndexedItem(
"layer",m_nGemHits,m_GemLayer);
765 status10 = m_nt10->addIndexedItem(
"mass",m_nGemHits,m_mass);
766 status10 = m_nt10->addIndexedItem(
"dchi2",m_nGemHits,m_Gchi2);
767 status10 = m_nt10->addIndexedItem(
"meas_phierr",m_nGemHits,m_meas_phierr);
768 status10 = m_nt10->addIndexedItem(
"meas_zerr",m_nGemHits,m_meas_zerr);
769 status10 = m_nt10->addIndexedItem(
"esti_phierr",m_nGemHits,m_esti_phierr);
770 status10 = m_nt10->addIndexedItem(
"esti_zerr",m_nGemHits,m_esti_zerr);
771 if( status10.isFailure() ) cout<<
"Ntuple10 add item failed!"<<endl;
777 NTuplePtr nt11(
ntupleSvc(),
"FILE104/n111");
779 if ( nt11 ) m_nt11 = nt11;
781 m_nt11=
ntupleSvc()->book(
"FILE104/n111",CLID_ColumnWiseTuple,
"truth");
783 status11 = m_nt11->addItem(
"evt",m_evt4);
784 status11 = m_nt11->addItem(
"ntruth",m_ntruth,0,400000);
785 status11 = m_nt11->addIndexedItem(
"phi",m_ntruth,m_dtphi);
786 status11 = m_nt11->addIndexedItem(
"z",m_ntruth,m_dtv);
787 status11 = m_nt11->addIndexedItem(
"postphi",m_ntruth,m_dtpostphi);
788 status11 = m_nt11->addIndexedItem(
"postz",m_ntruth,m_dtpostz);
789 status11 = m_nt11->addIndexedItem(
"layer",m_ntruth,m_tlayer);
790 if( status11.isFailure() ) cout<<
"Ntuple11 add item failed!"<<endl;
796 NTuplePtr nt12(
ntupleSvc(),
"FILE104/n112");
798 if ( nt12 ) m_nt12 = nt12;
800 m_nt12=
ntupleSvc()->book(
"FILE104/n112",CLID_ColumnWiseTuple,
"PatRecComp");
802 status12 = m_nt12->addItem(
"evt",m_evt);
803 status12 = m_nt12->addItem(
"track",m_track);
804 status12 = m_nt12->addItem(
"diff_dr",m_diff_dr);
805 status12 = m_nt12->addItem(
"diff_phi0",m_diff_phi0);
806 status12 = m_nt12->addItem(
"diff_kappa",m_diff_kappa);
807 status12 = m_nt12->addItem(
"diff_dz",m_diff_dz);
808 status12 = m_nt12->addItem(
"diff_tanl",m_diff_tanl);
809 status12 = m_nt12->addItem(
"diff_p",m_diff_p);
810 if( status12.isFailure() ) cout<<
"Ntuple12 add item failed!"<<endl;
923 MsgStream log(
msgSvc(), name());
924 log << MSG::INFO <<
"in execute()" << endreq;
925 if(
testOutput) std::cout<<
"begin to deal with EVENT ..."<<(++
eventno)<<std::endl;
971 StatusCode sc = service (
"MagneticFieldSvc",IMFSvc);
972 if(sc!=StatusCode::SUCCESS) {
973 log << MSG::ERROR <<
"Unable to open Magnetic field service"<<endreq;
982 std::cout<<
" execute, referField from MagneticFieldSvc: "<< (IMFSvc->
getReferField())*10000 <<std::endl;
989 IDataProviderSvc* evtSvc = NULL;
990 Gaudi::svcLocator()->service(
"EventDataSvc", evtSvc);
992 log << MSG::INFO <<
"makeTds:event Svc has been found" << endreq;
994 log << MSG::FATAL <<
"makeTds:Could not find eventSvc" << endreq;
995 return StatusCode::SUCCESS;
999 IDataManagerSvc *dataManSvc;
1000 dataManSvc=
dynamic_cast<IDataManagerSvc*
>(evtSvc);
1001 DataObject *aKalTrackCol;
1002 evtSvc->findObject(
"/Event/Recon/RecMdcKalTrackCol",aKalTrackCol);
1003 if(aKalTrackCol != NULL) {
1004 dataManSvc->clearSubTree(
"/Event/Recon/RecMdcKalTrackCol");
1005 evtSvc->unregisterObject(
"/Event/Recon/RecMdcKalTrackCol");
1008 kalsc = evtSvc->registerObject(
"/Event/Recon/RecMdcKalTrackCol", kalcol);
1009 if( kalsc.isFailure() ) {
1010 log << MSG::FATAL <<
"Could not register RecMdcKalTrack" << endreq;
1011 return StatusCode::SUCCESS;
1013 log << MSG::INFO <<
"RecMdcKalTrackCol registered successfully!" <<endreq;
1015 DataObject *aKalHelixSegCol;
1016 evtSvc->findObject(
"/Event/Recon/RecMdcKalHelixSegCol", aKalHelixSegCol);
1017 if(aKalHelixSegCol != NULL){
1018 dataManSvc->clearSubTree(
"/Event/Recon/RecMdcKalHelixSegCol");
1019 evtSvc->unregisterObject(
"/Event/Recon/RecMdcKalHelixSegCol");
1022 kalsc = evtSvc->registerObject(
"/Event/Recon/RecMdcKalHelixSegCol", helixsegcol);
1023 if( kalsc.isFailure()){
1024 log<< MSG::FATAL <<
"Could not register RecMdcKalHelixSeg" <<endreq;
1025 return StatusCode::SUCCESS;
1027 log << MSG::INFO <<
"RecMdcKalHelixSegCol register successfully!" <<endreq;
1037 ISvcLocator* svcLocator = Gaudi::svcLocator();
1039 StatusCode Cgem_sc=svcLocator->service(
"CgemGeomSvc", ISvc);
1041 if (!Cgem_sc.isSuccess()) log<< MSG::INFO <<
"KalFitAlg::execute(): Could not open CGEM geometry file" << endreq;
1047 std::cout<<
"ERROR OCCUR when dynamic_cast in KalFitAlg::execute ...!!"<<std::endl;
1050 SmartDataPtr<Event::EventHeader> eventHeader(eventSvc(),
"/Event/EventHeader");
1052 log << MSG::WARNING <<
"Could not find Event Header" << endreq;
1053 return StatusCode::FAILURE;
1055 int eventNo = eventHeader->eventNumber();
1057 int runNo = eventHeader->runNumber();
1061 if(
testOutput) cout<<endl<<
"$$$$$$$$$$$ run="<<
runNo<<
", evt="<<
eventNo<<
" $$$$$$$$$$$$$$$$$"<<endl<<endl;
1064 SmartDataPtr<RecEsTimeCol> estimeCol(eventSvc(),
"/Event/Recon/RecEsTimeCol");
1065 if (estimeCol && estimeCol->size()) {
1066 RecEsTimeCol::iterator iter_evt = estimeCol->begin();
1067 t0 = (*iter_evt)->getTest();
1070 log << MSG::WARNING <<
"Could not find EvTimeCol" << endreq;
1071 return StatusCode::SUCCESS;
1076 std::cout<<
"in KalFitAlg , we get the event start time = "<<t0<<std::endl;
1080 SmartDataPtr<MdcDigiCol> mdcDigiCol(evtSvc,
"/Event/Digi/MdcDigiCol");
1081 if (sc!=StatusCode::SUCCESS) {
1082 log << MSG::FATAL <<
"Could not find MdcDigiCol!" << endreq;
1083 return StatusCode::SUCCESS;
1094 m_evtid = eventHeader->eventNumber();
1097 SmartDataPtr<McParticleCol> mcPartCol(eventSvc(),
"/Event/MC/McParticleCol");
1099 log << MSG::WARNING <<
"Could not find McParticle" << endreq;
1104 McParticleCol::iterator i_mcTrk = mcPartCol->begin();
1105 for (;i_mcTrk != mcPartCol->end(); i_mcTrk++) {
1106 if(!(*i_mcTrk)->primaryParticle())
continue;
1107 const HepLorentzVector& mom((*i_mcTrk)->initialFourMomentum());
1108 const HepLorentzVector& pos = (*i_mcTrk)->initialPosition();
1109 log << MSG::DEBUG <<
"MCINFO:particleId=" << (*i_mcTrk)->particleProperty()
1110 <<
" theta=" << mom.theta() <<
" phi="<< mom.phi()
1111 <<
" px="<< mom.px() <<
" py="<< mom.py() <<
" pz="<< mom.pz()
1113 double charge = 0.0;
1114 int pid = (*i_mcTrk)->particleProperty();
1115 int mother_id = ((*i_mcTrk)->mother()).particleProperty();
1117 charge = m_particleTable->particle( pid )->charge();
1118 }
else if ( pid <0 ) {
1119 charge = m_particleTable->particle( -pid )->charge();
1122 log << MSG::WARNING <<
"wrong particle id, please check data" <<endreq;
1125 Hep3Vector mom2(mom.px(),mom.py(),mom.pz());
1128 log << MSG::DEBUG <<
"charge of the track " << charge << endreq;
1129 if(
debug_ == 4) cout<<
"helix: "<<mchelix.
a()<<endl;
1131 for(
int j =0; j<5; j++) {
1132 m_mchelix[j] = mchelix.
a()[j];
1135 m_mcptot = sqrt(1+pow(m_mchelix[4],2))/m_mchelix[2];
1138 cout<<
"MC pid, mother_id = "<<pid<<
", "<<mother_id<<endl;
1139 cout<<
" p4 = "<<(*i_mcTrk)->initialFourMomentum()<<endl;
1140 cout<<
" start from "<<(*i_mcTrk)->initialPosition()<<endl;
1141 cout<<
" end at "<<(*i_mcTrk)->finalPosition()<<endl;
1142 cout<<
" Helix: "<<mchelix.
a()<<endl;
1143 cout<<
"mc ptot, theta, phi, R = "<<m_mcptot<<
", "<<mom2.theta()/acos(-1.)*180<<
", "<<mom2.phi()/acos(-1.)*180<<
", "<<mchelix.
radius()<<endl;
1152 SmartDataPtr<RecMdcTrackCol> newtrkCol(eventSvc(),
"/Event/Recon/RecMdcTrackCol");
1154 log << MSG::FATAL <<
"Could not find RecMdcTrackCol" << endreq;
1155 return( StatusCode::SUCCESS);
1157 log << MSG::INFO <<
"Begin to make MdcRecTrkCol and MdcRecWirhitCol"<<endreq;
1163 mtrkadd_mgr->clear();
1167 double trkx1,trkx2,trky1,trky2,trkz1,trkz2,trkthe1,trkthe2,trkphi1,trkphi2,trkp1,trkp2,trkr1,trkr2,trkkap1,trkkap2,trktanl1,trktanl2;
1168 Hep3Vector csmp3[2];
1171 RecMdcTrackCol::iterator iter_trk = newtrkCol->begin();
1172 for(
int kj = 1; iter_trk != newtrkCol->end(); iter_trk++,kj++) {
1174 csmp3[kj-1]=(*iter_trk)->p3();
1175 csmphi[kj-1] = (*iter_trk)->phi();
1179 for(
int j = 0, ij = 0; j<5; j++) {
1180 m_trkhelix[j] = (*iter_trk)->helix()[j];
1182 for(
int k=0; k<=j; k++,ij++) {
1183 m_trkerror[ij] = (*iter_trk)->err()[j][k];
1187 m_trkptot = sqrt(1+pow(m_trkhelix[4],2))/m_trkhelix[2];
1189 m_trksigp = sqrt(pow((m_trkptot/m_trkhelix[2]),2)*m_trkerror[5]+
1190 pow((m_trkhelix[4]/m_trkptot),2)*pow((1/m_trkhelix[2]),4)*m_trkerror[14]-
1191 2*m_trkhelix[4]*m_trkerror[12]*pow((1/m_trkhelix[2]),3));
1193 m_trkndf = (*iter_trk)->ndof();
1194 m_trkchisq = (*iter_trk)->chi2();
1196 if (
debug_ == 4) cout<<
"Ea from RecMdcTrackCol..." <<(*iter_trk)->err()<<endl;
1198 StatusCode sc3 = m_nt3->write();
1199 if( sc3.isFailure() ) cout<<
"Ntuple3 filling failed!"<<endl;
1229 log << MSG::DEBUG <<
"retrieved MDC tracks:"
1230 <<
" Nhits " <<(*iter_trk)->getNhits()
1231 <<
" Nster " <<(*iter_trk)->nster() <<endreq;
1233 HitRefVec gothits = (*iter_trk)->getVecHits();
1237 rectrk->
id = (*iter_trk)->trackId();
1238 rectrk->
chiSq = (*iter_trk)->chi2();
1239 rectrk->
ndf = (*iter_trk)->ndof();
1240 rectrk->
fiTerm = (*iter_trk)->getFiTerm();
1241 rectrk->
nhits = (*iter_trk)->getNhits();
1242 rectrk->
nster = (*iter_trk)->nster();
1243 rectrk->
nclus = (*iter_trk)->getNcluster();
1244 rectrk->
stat = (*iter_trk)->stat();
1245 status_temp = (*iter_trk)->stat();
1247 trkadd->
id = (*iter_trk)->trackId();
1251 trkadd->
body = rectrk;
1252 rectrk->
add = trkadd;
1257 for (
int i=0; i<5; i++) {
1258 rectrk->
helix[i] = (*iter_trk)->helix()[i];
1259 if( i<3 ) rectrk->
pivot[i] = (*iter_trk)->getPivot()[i];
1260 for(
int j = 1; j<i+2;j++) {
1261 rectrk->
error[i*(i+1)/2+j-1] = (*iter_trk)->err()(i+1,j);
1264 std::sort(gothits.begin(), gothits.end(), order_rechits);
1265 HitRefVec::iterator it_gothit = gothits.begin();
1266 for( ; it_gothit != gothits.end(); it_gothit++) {
1268 if( (*it_gothit)->getStat() != 1 ) {
1270 log<<MSG::WARNING<<
"this hit is not used in helix fitting!"<<endreq;
1275 log << MSG::DEBUG <<
"retrieved hits in MDC tracks:"
1276 <<
" hits DDL " <<(*it_gothit)->getDriftDistLeft()
1277 <<
" hits DDR " <<(*it_gothit)->getDriftDistRight()
1278 <<
" error DDL " <<(*it_gothit)->getErrDriftDistLeft()
1279 <<
" error DDR " <<(*it_gothit)->getErrDriftDistRight()
1280 <<
" id of hit "<<(*it_gothit)->getId()
1281 <<
" track id of hit "<<(*it_gothit)->getTrkId()
1282 <<
" hits ADC " <<(*it_gothit)->getAdc() << endreq;
1285 whit->
id = (*it_gothit)->getId();
1286 whit->
ddl = (*it_gothit)->getDriftDistLeft();
1287 whit->
ddr = (*it_gothit)->getDriftDistRight();
1288 whit->
erddl = (*it_gothit)->getErrDriftDistLeft();
1289 whit->
erddr = (*it_gothit)->getErrDriftDistRight();
1290 whit->
pChiSq = (*it_gothit)->getChisqAdd();
1291 whit->
lr = (*it_gothit)->getFlagLR();
1292 whit->
stat = (*it_gothit)->getStat();
1293 mdcid = (*it_gothit)->getMdcId();
1297 int wid = w0id + localwid;
1299 <<
"lr from PR: "<<whit->
lr
1300 <<
" layerId = " << layid
1301 <<
" wireId = " << localwid
1311 whit->
tdc = (*it_gothit)->getTdc();
1312 whit->
adc= (*it_gothit)->getAdc();
1313 rectrk->
hitcol.push_back(whit);
1314 mhit_mgr->push_back(*whit);
1316 mtrk_mgr->push_back(*rectrk);
1317 mtrkadd_mgr->push_back(*trkadd);
1325 m_trkdelx = trkx1 - trkx2;
1326 m_trkdely = trky1 - trky2;
1327 m_trkdelz = trkz1 - trkz2;
1328 m_trkdelthe = trkthe1 + trkthe2;
1329 m_trkdelphi = trkphi1- trkphi2;
1330 m_trkdelp = trkp1 - trkp2;
1331 StatusCode sc4 = m_nt4->write();
1332 if( sc4.isFailure() ) cout<<
"Ntuple4 filling failed!"<<endl;
1335 if(
debug_ == 4) { std::cout<<
"before refit,ntrk,nhits,nadd..."<<mtrk_mgr->size()
1336 <<
"********"<<mhit_mgr->size()<<
"****"<<mtrkadd_mgr->size()<<endl;
1347 double mdang = 180.0 - csmp3[0].angle(csmp3[1].
unit())*180.0/
M_PI;
1348 double mdphi = 180.0 - fabs(csmphi[0]-csmphi[1])*180.0/
M_PI;
1353 log << MSG::DEBUG <<
"after kalman_fitting(),but in execute...."<<endreq;
1360 SmartDataPtr<RecMdcKalTrackCol> recmdckaltrkCol(eventSvc(),
"/Event/Recon/RecMdcKalTrackCol");
1365 RecMdcKalTrackCol::iterator KalTrk= recmdckaltrkCol->begin();
1367 for(;KalTrk !=recmdckaltrkCol->end();KalTrk++, i_trk++){
1369 bool trkFailed =
false;
1370 for(
int hypo=0; hypo<5; hypo++)
1372 if((*KalTrk)->getStat(0,hypo)==1) {
1373 nFailedTrks[hypo]++;
1377 if(
debug_==5&&trkFailed) cout<<
" Evt "<<
myEventNo<<
", track "<<i_trk<<
": Kalman filter failed"<<endl;
1380 HelixSegRefVec::iterator iter_hit = gothelixsegs.begin();
1382 int nhitofthistrk=0;
1383 for( ; iter_hit != gothelixsegs.end(); iter_hit++){
1390 iter_hit=gothelixsegs.begin();
1438 SmartDataPtr<Event::CgemMcHitCol> cgemMcHitCol(eventSvc(),
"/Event/MC/CgemMcHitCol");
1440 log << MSG::WARNING <<
"Could not retrieve Cgem MC truth" << endreq;
1441 return StatusCode::FAILURE;
1443 m_evt4=eventHeader->eventNumber();
1444 Event::CgemMcHitCol::iterator iter_truth = cgemMcHitCol->begin();
1446 double dprex,dprey,dprez,dpostx,dposty,dpostz;
1448 for(; iter_truth != cgemMcHitCol->end(); ++iter_truth){
1449 layer = (*iter_truth)->GetLayerID();
1450 dprex = (*iter_truth)->GetPositionXOfPrePoint();
1451 dprey = (*iter_truth)->GetPositionYOfPrePoint();
1452 dprez = (*iter_truth)->GetPositionZOfPrePoint();
1453 dpostx = (*iter_truth)->GetPositionXOfPostPoint();
1454 dposty = (*iter_truth)->GetPositionYOfPostPoint();
1455 dpostz = (*iter_truth)->GetPositionZOfPostPoint();
1463 if((dprex >=0&&dprey>=0)||(dprex <0&&dprey>=0)){
1465 m_dtphi[jj] = acos(dprex/diR);
1468 m_dtpostphi[jj] = acos(dpostx/doR);
1469 m_dtpostz[jj]= dpostz;
1471 if((dprex <0&&dprey<0)||(dprex >=0&&dprey<0)){
1473 m_dtphi[jj] = -acos(dprex/diR);
1476 m_dtpostphi[jj] = -acos(dpostx/doR);
1477 m_dtpostz[jj]= dpostz;
1479 midx=(dprex+dpostx)/2;
1480 midy=(dprey+dposty)/2;
1481 if((midx>=0&&midy>=0)||(midx<0&&midy>=0)){
1483 m_dtphi[jj]=acos(midx/dmR);
1485 if((midx<0&&midy<0)||(midx>=0&&midy<0)){
1487 m_dtphi[jj]=-acos(midx/dmR);
1489 m_dtv[jj] = (m_dtv[jj]+m_dtpostz[jj])/20;
1490 m_tlayer[jj] = layer;
1500 return StatusCode::SUCCESS;
3198 cout<<
"**********************"<<endl;
3199 cout<<
"filter pid type "<<l_mass<<endl;
3204 StatusCode sc = Gaudi::svcLocator()->service(
"MdcGeomSvc", igeomsvc);
3205 if(sc==StatusCode::FAILURE) cout <<
"GeoSvc failing!!!!!!!SC=" << sc << endl;
3208 std::cout<<
"ERROR OCCUR when dynamic_cast in KalFitTrack.cxx !!"<<std::endl;
3212 HepVector pos_old(3,0);
3213 double r0kal_prec(0);
3215 int nhit = track.
HitsMdc().size();
3216 if(
debug_ == 4) cout<<
"filter_fwd..........111 nhit="<<nhit<<endl;
3217 for(
int i=0 ; i < nhit; i++ ) {
3220 if (HitMdc.
chi2()<0)
continue;
3226 int wireid = Wire.
geoID();
3230 Hep3Vector wire = (CLHEP::Hep3Vector)fwd -(CLHEP::Hep3Vector)bck;
3233 work.
pivot((fwd + bck) * .5);
3240 HepPoint3D x0kal = (work.
x(0).z() - bck.z())/ wire.z() * wire + bck;
3243 std::cout<<
" a="<<track.
a()<<std::endl;
3244 std::cout<<
" pivot="<<track.
pivot()<<std::endl;
3245 std::cout<<
" x0kal before sag: "<<x0kal<<std::endl;
3290 double tension = geowire->
Tension();
3292 double zinit(x0kal.z()), lzx(Wire.
lzx());
3294 double A = 47.35E-6/tension;
3295 double Zp = (zinit - bck.z())*lzx/wire.z();
3298 std::cout<<
" sag in filter_fwd_anal: "<<std::endl;
3299 std::cout<<
" x0kal.x(): "<<std::setprecision(10)<<x0kal.x()<<std::endl;
3300 std::cout<<
"zinit: "<<zinit<<
" bck.z(): "<<bck.z()<<std::endl;
3301 std::cout<<
" wire.x()*(zinit-bck.z())/wire.z(): "<<std::setprecision(10)
3302 <<(wire.x()*(zinit-bck.z())/wire.z())<<std::endl;
3303 std::cout<<
"bck.x(): "<<std::setprecision(10)<<bck.x()<<std::endl;
3304 std::cout<<
" wire.x()*(zinit-bck.z())/wire.z() + bck.x(): "<<std::setprecision(10)
3305 <<(wire.x()*(zinit-bck.z())/wire.z() + bck.x())<<std::endl;
3308 result.setX(wire.x()*(zinit-bck.z())/wire.z() + bck.x());
3309 result.setY((A*(Zp-lzx)+wire.y()/lzx)*Zp+bck.y());
3310 result.setZ((A*(2*Zp-lzx)*lzx+wire.y())/wire.z());
3312 wire.setX(wire.x()/wire.z());
3313 wire.setY(result.z());
3315 x0kal.setX(result.x());
3316 x0kal.setY(result.y());
3321 double r0kal = x0kal.perp();
3326 if(i==0) track.
pivot(x0kal);
3330 if (nhits_read == 1) {
3341 std::cout<<
" x0kal after sag: "<<x0kal<<std::endl;
3342 std::cout<<
" a="<<track.
a()<<std::endl;
3343 std::cout<<
" Ea="<<track.
Ea()<<std::endl;
3346 double dtracknew = 0.;
3350 if(fabs(track.
kappa())>0&&fabs(track.
kappa())<1000.0&&fabs(track.
tanl())<7.02) {
3351 Hep3Vector meas = track.
momentum(0).cross(wire).unit();
3352 double diff_chi2 = track.
chiSq();
3353 Hep3Vector IP(0,0,0);
3359 for(
unsigned k=i+1 ; k < nhit; k++ )
3364 double dchi2 = -1.0;
3419 m_residest = dtrack-dtdc;
3420 m_residnew = dtracknew -dtdc;
3423 m_anal_dr = worktemp.
a()[0];
3424 m_anal_phi0 = worktemp.
a()[1];
3425 m_anal_kappa = worktemp.
a()[2];
3426 m_anal_dz = worktemp.
a()[3];
3427 m_anal_tanl = worktemp.
a()[4];
3428 m_anal_ea_dr = worktemp.
Ea()[0][0];
3429 m_anal_ea_phi0 = worktemp.
Ea()[1][1];
3430 m_anal_ea_kappa = worktemp.
Ea()[2][2];
3431 m_anal_ea_dz = worktemp.
Ea()[3][3];
3432 m_anal_ea_tanl = worktemp.
Ea()[4][4];
3433 StatusCode sc5 = m_nt5->write();
3434 if( sc5.isFailure() ) cout<<
"Ntuple5 filling failed!"<<endl;
3442 m_dchi2_hit_e[m_hit_no] = dchi2;
3443 m_residest_hit_e[m_hit_no] = dtrack-dtdc;
3444 m_residnew_hit_e[m_hit_no] = dtracknew -dtdc;
3447 m_anal_dr_hit_e[m_hit_no] = worktemp.
a()[0];
3448 m_anal_phi0_hit_e[m_hit_no] = worktemp.
a()[1];
3449 m_anal_kappa_hit_e[m_hit_no] = worktemp.
a()[2];
3450 m_anal_dz_hit_e[m_hit_no] = worktemp.
a()[3];
3451 m_anal_tanl_hit_e[m_hit_no] = worktemp.
a()[4];
3452 m_anal_ea_dr_hit_e[m_hit_no] = worktemp.
Ea()[0][0];
3453 m_anal_ea_phi0_hit_e[m_hit_no] = worktemp.
Ea()[1][1];
3454 m_anal_ea_kappa_hit_e[m_hit_no] = worktemp.
Ea()[2][2];
3455 m_anal_ea_dz_hit_e[m_hit_no] = worktemp.
Ea()[3][3];
3456 m_anal_ea_tanl_hit_e[m_hit_no] = worktemp.
Ea()[4][4];
3460 else if(l_mass == 1){
3461 m_dchi2_hit_mu[m_hit_no] = dchi2;
3462 m_residest_hit_mu[m_hit_no] = dtrack-dtdc;
3463 m_residnew_hit_mu[m_hit_no] = dtracknew -dtdc;
3466 m_anal_dr_hit_mu[m_hit_no] = worktemp.
a()[0];
3467 m_anal_phi0_hit_mu[m_hit_no] = worktemp.
a()[1];
3468 m_anal_kappa_hit_mu[m_hit_no] = worktemp.
a()[2];
3469 m_anal_dz_hit_mu[m_hit_no] = worktemp.
a()[3];
3470 m_anal_tanl_hit_mu[m_hit_no] = worktemp.
a()[4];
3471 m_anal_ea_dr_hit_mu[m_hit_no] = worktemp.
Ea()[0][0];
3472 m_anal_ea_phi0_hit_mu[m_hit_no] = worktemp.
Ea()[1][1];
3473 m_anal_ea_kappa_hit_mu[m_hit_no] = worktemp.
Ea()[2][2];
3474 m_anal_ea_dz_hit_mu[m_hit_no] = worktemp.
Ea()[3][3];
3475 m_anal_ea_tanl_hit_mu[m_hit_no] = worktemp.
Ea()[4][4];
3480 else if(l_mass == 2){
3481 m_dchi2_hit_pi[m_hit_no] = dchi2;
3482 m_residest_hit_pi[m_hit_no] = dtrack-dtdc;
3483 m_residnew_hit_pi[m_hit_no] = dtracknew -dtdc;
3486 m_anal_dr_hit_pi[m_hit_no] = worktemp.
a()[0];
3487 m_anal_phi0_hit_pi[m_hit_no] = worktemp.
a()[1];
3488 m_anal_kappa_hit_pi[m_hit_no] = worktemp.
a()[2];
3489 m_anal_dz_hit_pi[m_hit_no] = worktemp.
a()[3];
3490 m_anal_tanl_hit_pi[m_hit_no] = worktemp.
a()[4];
3491 m_anal_ea_dr_hit_pi[m_hit_no] = worktemp.
Ea()[0][0];
3492 m_anal_ea_phi0_hit_pi[m_hit_no] = worktemp.
Ea()[1][1];
3493 m_anal_ea_kappa_hit_pi[m_hit_no] = worktemp.
Ea()[2][2];
3494 m_anal_ea_dz_hit_pi[m_hit_no] = worktemp.
Ea()[3][3];
3495 m_anal_ea_tanl_hit_pi[m_hit_no] = worktemp.
Ea()[4][4];
3500 else if(l_mass == 3){
3501 m_dchi2_hit_k[m_hit_no] = dchi2;
3502 m_residest_hit_k[m_hit_no] = dtrack-dtdc;
3503 m_residnew_hit_k[m_hit_no] = dtracknew -dtdc;
3506 m_anal_dr_hit_k[m_hit_no] = worktemp.
a()[0];
3507 m_anal_phi0_hit_k[m_hit_no] = worktemp.
a()[1];
3508 m_anal_kappa_hit_k[m_hit_no] = worktemp.
a()[2];
3509 m_anal_dz_hit_k[m_hit_no] = worktemp.
a()[3];
3510 m_anal_tanl_hit_k[m_hit_no] = worktemp.
a()[4];
3511 m_anal_ea_dr_hit_k[m_hit_no] = worktemp.
Ea()[0][0];
3512 m_anal_ea_phi0_hit_k[m_hit_no] = worktemp.
Ea()[1][1];
3513 m_anal_ea_kappa_hit_k[m_hit_no] = worktemp.
Ea()[2][2];
3514 m_anal_ea_dz_hit_k[m_hit_no] = worktemp.
Ea()[3][3];
3515 m_anal_ea_tanl_hit_k[m_hit_no] = worktemp.
Ea()[4][4];
3520 else if(l_mass == 4){
3521 m_dchi2_hit_p[m_hit_no] = dchi2;
3522 m_residest_hit_p[m_hit_no] = dtrack-dtdc;
3523 m_residnew_hit_p[m_hit_no] = dtracknew -dtdc;
3526 m_anal_dr_hit_p[m_hit_no] = worktemp.
a()[0];
3527 m_anal_phi0_hit_p[m_hit_no] = worktemp.
a()[1];
3528 m_anal_kappa_hit_p[m_hit_no] = worktemp.
a()[2];
3529 m_anal_dz_hit_p[m_hit_no] = worktemp.
a()[3];
3530 m_anal_tanl_hit_p[m_hit_no] = worktemp.
a()[4];
3531 m_anal_ea_dr_hit_p[m_hit_no] = worktemp.
Ea()[0][0];
3532 m_anal_ea_phi0_hit_p[m_hit_no] = worktemp.
Ea()[1][1];
3533 m_anal_ea_kappa_hit_p[m_hit_no] = worktemp.
Ea()[2][2];
3534 m_anal_ea_dz_hit_p[m_hit_no] = worktemp.
Ea()[3][3];
3535 m_anal_ea_tanl_hit_p[m_hit_no] = worktemp.
Ea()[4][4];
3546 diff_chi2 = chi2 - diff_chi2;
3547 HitMdc.
chi2(diff_chi2);
3857 MsgStream log(
msgSvc(), name());
3858 double Pt_threshold(0.3);
3859 Hep3Vector IP(0,0,0);
3865 if ( !&whMgr )
return;
3868 int ntrk = mdcMgr->size();
3869 double* rPt =
new double[ntrk];
3870 int* rOM =
new int[ntrk];
3871 unsigned int* order =
new unsigned int[ntrk];
3872 unsigned int* rCont =
new unsigned int[ntrk];
3873 unsigned int* rGen =
new unsigned int[ntrk];
3875 if(
testOutput) cout<<
"nMdcTrk = "<<ntrk<<endl;
3878 for(vector<MdcRec_trk>::iterator it = mdcMgr->begin(),
3879 end = mdcMgr->end(); it != end; it++) {
3883 rPt[index] = 1 / fabs(it->helix[2]);
3884 if(
debug_ == 4) cout<<
"rPt...[ "<<index<<
" ]...."<< rPt[index] <<endl;
3885 if(rPt[index] < 0) rPt[index] =
DBL_MAX;
3887 std::vector<MdcRec_wirhit*> pt = it->hitcol ;
3888 if(
debug_ == 4) cout<<
"ppt size: "<< pt.size()<<endl;
3890 for (vector<MdcRec_wirhit*>::iterator ii = pt.end()-1;
3891 ii !=pt.begin()-1; ii--) {
3892 int lyr((*ii)->geo->Lyr()->Id());
3893 if (outermost < lyr) outermost = lyr;
3894 if(
debug_ == 4) cout<<
"outmost: "<<outermost<<
" lyr: "<<lyr<<endl;
3896 rOM[index] = outermost;
3897 order[index] = index;
3902 for (
int j, k = ntrk - 1; k >= 0; k = j){
3904 for(
int i = 1; i <= k; i++)
3905 if(rPt[i - 1] < rPt[i]){
3907 std::swap(order[i], order[j]);
3908 std::swap(rPt[i], rPt[j]);
3909 std::swap(rOM[i], rOM[j]);
3910 std::swap(rCont[i], rCont[j]);
3911 std::swap(rGen[i], rGen[j]);
3919 DataObject *aReconEvent;
3920 eventSvc()->findObject(
"/Event/Recon",aReconEvent);
3924 StatusCode sc = eventSvc()->registerObject(
"/Event/Recon",recevt );
3925 if(sc!=StatusCode::SUCCESS) {
3926 log << MSG::FATAL <<
"Could not register ReconEvent" <<endreq;
3934 log << MSG::INFO <<
"beginning to make RecMdcKalTrackCol" <<endreq;
3936 for(
int l = 0; l < ntrk; l++) {
3945 MdcRec_trk& TrasanTRK = *(mdcMgr->begin() + order[l]);
3949 int trasqual = TrasanTRK_add.
quality;
3950 if(
debug_ == 4) cout<<
"kalman_fitting>trasqual... "<<trasqual<<endl;
3951 if (trasqual)
continue;
3955 cout <<
"******* KalFit NUMBER : " << newcount << std::endl;
3959 if ((TrasanTRK_add.
decision & 32) == 32 ||
3960 (TrasanTRK_add.
decision & 64) == 64) type = 1;
3965 TrasanTRK.
pivot[2]);
3968 for(
int i = 0; i < 5; i++)
3969 a[i] = TrasanTRK.
helix[i];
3972 for(
int i = 0, k = 0; i < 5; i++) {
3973 for(
int j = 0; j <= i; j++) {
3975 ea[j][i] = ea[i][j];
3978 double fiTerm = TrasanTRK.
fiTerm;
3984 double pT_trk = fabs(track_lead.
pt());
3987 int inlyr(999), outlyr(-1);
3988 int* rStat =
new int[43];
3989 for(
int irStat=0;irStat<43;++irStat)rStat[irStat]=0;
3990 std::vector<MdcRec_wirhit*> pt=TrasanTRK.
hitcol;
3992 if(
debug_ == 4) cout<<
"*********Pt size****"<< pt.size()<<endl;
3995 for (vector<MdcRec_wirhit*>::iterator ii = pt.end()-1;
3996 ii != pt.begin()-1; ii--) {
3997 Num[(*ii)->geo->Lyr()->Id()]++;
4000 for (vector<MdcRec_wirhit*>::iterator ii = pt.end()-1;
4001 ii != pt.begin()-1; ii--) {
4005 if( (useNCGem_==0&&Num[(*ii)->geo->Lyr()->Id()]>3) || (useNCGem_>0&&pT_trk>0.12 && Num[(*ii)->geo->Lyr()->Id()]>3) )
4008 cout <<
"WARNING: I found " << Num[(*ii)->geo->Lyr()->Id()]
4009 <<
" hits in the layer "
4010 << (*ii)->geo->Lyr()->Id() << std::endl;
4030 double dist[2] = {rechit.
ddl, rechit.
ddr};
4031 double erdist[2] = {rechit.
erddl, rechit.
erddr};
4036 if (rechit.
lr==2 || rechit.
lr==0) lr_decision=-1;
4038 else if (rechit.
lr==1) lr_decision=1;
4041 int ind(geo->
Lyr()->
Id());
4045 lr_decision, rechit.
tdc,
4050 if (inlyr>ind) inlyr = ind;
4051 if (outlyr<ind) outlyr = ind;
4054 cout <<
"**** NUMBER OF Mdc HITS (TRASAN) = " << hit_asso << std::endl;
4057 int empty_between(0), empty(0);
4058 for (
int i= inlyr; i <= outlyr; i++)
4059 if (!rStat[i]) empty_between++;
4060 empty = empty_between+inlyr+(42-outlyr);
4067 for(std::vector<KalFitHitMdc>::iterator it_hit = track_lead.
HitsMdc().begin(); it_hit!=track_lead.
HitsMdc().end(); it_hit++){
4073 track_lead.
type(type);
4075 unsigned int nhit = track_lead.
HitsMdc().size();
4076 if (!nhit &&
debug_ == 4) {
4077 cout <<
" ATTENTION TRACK WITH ONLY HITS " << nhit << std::endl;
4082 double KalFitst(0), KalFitax(0), KalFitschi2(0);
4085 Hep3Vector outer_pivot(track_lead.
x(fiTerm));
4086 track_lead.
pivot(outer_pivot);
4092 HepSymMatrix Eakal(5,0);
4095 double costheta = track_lead.
a()[4] / sqrt(1.0 + track_lead.
a()[4]*track_lead.
a()[4]);
4105 cout <<
"from Mdc Pattern Recognition: " << std::endl;
4111 cout <<
" dr = " << work.
a()[0]
4112 <<
", Er_dr = " << sqrt(work.
Ea()[0][0]) << std::endl;
4113 cout <<
" phi0 = " << work.
a()[1]
4114 <<
", Er_phi0 = " << sqrt(work.
Ea()[1][1]) << std::endl;
4115 cout <<
" PT = " << 1/work.
a()[2]
4116 <<
", Er_kappa = " << sqrt(work.
Ea()[2][2]) << std::endl;
4117 cout <<
" dz = " << work.
a()[3]
4118 <<
", Er_dz = " << sqrt(work.
Ea()[3][3]) << std::endl;
4119 cout <<
" tanl = " << work.
a()[4]
4120 <<
", Er_tanl = " << sqrt(work.
Ea()[4][4]) << std::endl;
4124 cout<<
"----------------------"<<endl;
4125 cout<<
"track "<<l<<
", pid = "<<
lead_<<
": "<<endl;
4142 cout <<
" Mdc FIRST KALMAN FIT " << std::endl;
4147 cout <<
" dr = " << work.
a()[0]
4148 <<
", Er_dr = " << sqrt(work.
Ea()[0][0]) << std::endl;
4149 cout <<
" phi0 = " << work.
a()[1]
4150 <<
", Er_phi0 = " << sqrt(work.
Ea()[1][1]) << std::endl;
4151 cout <<
" PT = " << 1/work.
a()[2]
4152 <<
", Er_kappa = " << sqrt(work.
Ea()[2][2]) << std::endl;
4153 cout <<
" dz = " << work.
a()[3]
4154 <<
", Er_dz = " << sqrt(work.
Ea()[3][3]) << std::endl;
4155 cout <<
" tanl = " << work.
a()[4]
4156 <<
", Er_tanl = " << sqrt(work.
Ea()[4][4]) << std::endl;
4163 complete_track(TrasanTRK, TrasanTRK_add, track_lead, kaltrk,kalcol,segcol,1);
4167 IDataProviderSvc* eventSvc = NULL;
4168 Gaudi::svcLocator()->service(
"EventDataSvc", eventSvc);
4170 log << MSG::INFO <<
"makeTds:event Svc has been found" << endreq;
4172 log << MSG::FATAL <<
"makeTds:Could not find eventSvc" << endreq;
4176 IDataManagerSvc *dataManSvc;
4177 dataManSvc=
dynamic_cast<IDataManagerSvc*
>(eventSvc);
4178 DataObject *aKalTrackCol;
4179 eventSvc->findObject(
"/Event/Recon/RecMdcKalTrackCol",aKalTrackCol);
4180 if(aKalTrackCol != NULL) {
4181 dataManSvc->clearSubTree(
"/Event/Recon/RecMdcKalTrackCol");
4182 eventSvc->unregisterObject(
"/Event/Recon/RecMdcKalTrackCol");
4184 kalsc = eventSvc->registerObject(
"/Event/Recon/RecMdcKalTrackCol", kalcol);
4185 if( kalsc.isFailure() ) {
4186 log << MSG::FATAL <<
"Could not register RecMdcKalTrack" << endreq;
4189 log << MSG::INFO <<
"RecMdcKalTrackCol registered successfully!" <<endreq;
4192 DataObject *aRecKalSegEvent;
4193 eventSvc->findObject(
"/Event/Recon/RecMdcKalHelixSegCol", aRecKalSegEvent);
4194 if(aRecKalSegEvent!=NULL) {
4196 segsc = eventSvc->unregisterObject(
"/Event/Recon/RecMdcKalHelixSegCol");
4197 if(segsc != StatusCode::SUCCESS) {
4198 log << MSG::FATAL <<
"Could not unregister RecMdcKalHelixSegCol collection" << endreq;
4202 segsc = eventSvc->registerObject(
"/Event/Recon/RecMdcKalHelixSegCol", segcol);
4203 if( segsc.isFailure() ) {
4204 log << MSG::FATAL <<
"Could not register RecMdcKalHelixSeg" << endreq;
4207 log << MSG::INFO <<
"RecMdcKalHelixSegCol registered successfully!" <<endreq;
4209 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.);
4210 double r1(0.),r2(0.),kap1(999.),kap2(999.),tanl1(0),tanl2(0.);
4211 double px1(0.),px2(0.),py1(0.),py2(0.),pz1(0.),pz2(0.),charge1(0.),charge2(0.);
4214 SmartDataPtr<RecMdcKalTrackCol> kaltrkCol(eventSvc,
"/Event/Recon/RecMdcKalTrackCol");
4216 log << MSG::FATAL <<
"Could not find RecMdcKalTrackCol" << endreq;
4219 log << MSG::INFO <<
"Begin to check RecMdcKalTrackCol"<<endreq;
4220 RecMdcKalTrackCol::iterator iter_trk = kaltrkCol->begin();
4221 for(
int jj=1; iter_trk != kaltrkCol->end(); iter_trk++,jj++) {
4222 log << MSG::DEBUG <<
"retrieved MDC Kalmantrack:"
4223 <<
"Track Id: " << (*iter_trk)->getTrackId()
4224 <<
" Mass of the fit: "<< (*iter_trk)->getMass(2)<< endreq
4225 <<
" Length of the track: "<< (*iter_trk)->getLength(2)
4226 <<
" Tof of the track: "<< (*iter_trk)->getTof(2) << endreq
4227 <<
" Chisq of the fit: "<< (*iter_trk)->getChisq(0,2)
4228 <<
" "<< (*iter_trk)->getChisq(1,2) << endreq
4229 <<
"Ndf of the fit: "<< (*iter_trk)->getNdf(0,2)
4230 <<
" "<< (*iter_trk)->getNdf(1,2) << endreq
4231 <<
"Kappa " << (*iter_trk)->getZHelix()[2]
4233 for(
int i = 0; i<43; i++) {
4234 log << MSG::DEBUG <<
"retrieved pathl["<<i<<
"]= "
4235 << (*iter_trk)->getPathl(i) <<endreq;
4238 m_trackid = (*iter_trk)->getTrackId();
4240 for(
int jj =0, iii=0; jj<5; jj++) {
4242 m_length[jj] = (*iter_trk)->getLength(jj);
4243 m_tof[jj] = (*iter_trk)->getTof(jj);
4244 m_nhits[jj] = (*iter_trk)->getNhits(jj);
4245 m_zhelix[jj] = (*iter_trk)->getZHelix()[jj];
4246 m_zhelixe[jj] = (*iter_trk)->getZHelixE()[jj];
4247 m_zhelixmu[jj] = (*iter_trk)->getZHelixMu()[jj];
4248 m_zhelixk[jj] = (*iter_trk)->getZHelixK()[jj];
4249 m_zhelixp[jj] = (*iter_trk)->getZHelixP()[jj];
4250 m_fhelix[jj] = (*iter_trk)->getFHelix()[jj];
4251 m_fhelixe[jj] = (*iter_trk)->getFHelixE()[jj];
4252 m_fhelixmu[jj] = (*iter_trk)->getFHelixMu()[jj];
4253 m_fhelixk[jj] = (*iter_trk)->getFHelixK()[jj];
4254 m_fhelixp[jj] = (*iter_trk)->getFHelixP()[jj];
4255 m_lhelix[jj] = (*iter_trk)->getLHelix()[jj];
4256 m_lhelixe[jj] = (*iter_trk)->getLHelixE()[jj];
4257 m_lhelixmu[jj] = (*iter_trk)->getLHelixMu()[jj];
4258 m_lhelixk[jj] = (*iter_trk)->getLHelixK()[jj];
4259 m_lhelixp[jj] = (*iter_trk)->getLHelixP()[jj];
4262 for(
int kk=0; kk<=jj; kk++,iii++) {
4263 m_zerror[iii] = (*iter_trk)->getZError()[jj][kk];
4264 m_zerrore[iii] = (*iter_trk)->getZErrorE()[jj][kk];
4265 m_zerrormu[iii] = (*iter_trk)->getZErrorMu()[jj][kk];
4266 m_zerrork[iii] = (*iter_trk)->getZErrorK()[jj][kk];
4267 m_zerrorp[iii] = (*iter_trk)->getZErrorP()[jj][kk];
4268 m_ferror[iii] = (*iter_trk)->getFError()[jj][kk];
4269 m_ferrore[iii] = (*iter_trk)->getFErrorE()[jj][kk];
4270 m_ferrormu[iii] = (*iter_trk)->getFErrorMu()[jj][kk];
4271 m_ferrork[iii] = (*iter_trk)->getFErrorK()[jj][kk];
4272 m_ferrorp[iii] = (*iter_trk)->getFErrorP()[jj][kk];
4273 m_lerror[iii] = (*iter_trk)->getLError()[jj][kk];
4274 m_lerrore[iii] = (*iter_trk)->getLErrorE()[jj][kk];
4275 m_lerrormu[iii] = (*iter_trk)->getLErrorMu()[jj][kk];
4276 m_lerrork[iii] = (*iter_trk)->getLErrorK()[jj][kk];
4277 m_lerrorp[iii] = (*iter_trk)->getLErrorP()[jj][kk];
4317 m_chisq[0][0] = (*iter_trk)->getChisq(0,0);
4318 m_chisq[1][0] = (*iter_trk)->getChisq(0,1);
4319 m_chisq[2][0] = (*iter_trk)->getChisq(0,2);
4320 m_chisq[3][0] = (*iter_trk)->getChisq(0,3);
4321 m_chisq[4][0] = (*iter_trk)->getChisq(0,4);
4322 m_chisq[0][1] = (*iter_trk)->getChisq(1,0);
4323 m_chisq[1][1] = (*iter_trk)->getChisq(1,1);
4324 m_chisq[2][1] = (*iter_trk)->getChisq(1,2);
4325 m_chisq[3][1] = (*iter_trk)->getChisq(1,3);
4326 m_chisq[4][1] = (*iter_trk)->getChisq(1,4);
4328 m_ndf[0][0] = (*iter_trk)->getNdf(0,0);
4329 m_ndf[1][0] = (*iter_trk)->getNdf(0,1);
4330 m_ndf[2][0] = (*iter_trk)->getNdf(0,2);
4331 m_ndf[3][0] = (*iter_trk)->getNdf(0,3);
4332 m_ndf[4][0] = (*iter_trk)->getNdf(0,4);
4333 m_ndf[0][1] = (*iter_trk)->getNdf(1,0);
4334 m_ndf[1][1] = (*iter_trk)->getNdf(1,1);
4335 m_ndf[2][1] = (*iter_trk)->getNdf(1,2);
4336 m_ndf[3][1] = (*iter_trk)->getNdf(1,3);
4337 m_ndf[4][1] = (*iter_trk)->getNdf(1,4);
4339 m_stat[0][0] = (*iter_trk)->getStat(0,0);
4340 m_stat[1][0] = (*iter_trk)->getStat(0,1);
4341 m_stat[2][0] = (*iter_trk)->getStat(0,2);
4342 m_stat[3][0] = (*iter_trk)->getStat(0,3);
4343 m_stat[4][0] = (*iter_trk)->getStat(0,4);
4344 m_stat[0][1] = (*iter_trk)->getStat(1,0);
4345 m_stat[1][1] = (*iter_trk)->getStat(1,1);
4346 m_stat[2][1] = (*iter_trk)->getStat(1,2);
4347 m_stat[3][1] = (*iter_trk)->getStat(1,3);
4348 m_stat[4][1] = (*iter_trk)->getStat(1,4);
4350 m_fptot = sqrt(1+pow(m_fhelix[4],2))/m_fhelix[2];
4351 m_fptote = sqrt(1+pow(m_fhelixe[4],2))/m_fhelixe[2];
4352 m_fptotmu = sqrt(1+pow(m_fhelixmu[4],2))/m_fhelixmu[2];
4353 m_fptotk = sqrt(1+pow(m_fhelixk[4],2))/m_fhelixk[2];
4354 m_fptotp = sqrt(1+pow(m_fhelixp[4],2))/m_fhelixp[2];
4356 m_lptot = sqrt(1+pow(m_lhelix[4],2))/m_lhelix[2];
4357 m_lptote = sqrt(1+pow(m_lhelixe[4],2))/m_lhelixe[2];
4358 m_lptotmu = sqrt(1+pow(m_lhelixmu[4],2))/m_lhelixmu[2];
4359 m_lptotk = sqrt(1+pow(m_lhelixk[4],2))/m_lhelixk[2];
4360 m_lptotp = sqrt(1+pow(m_lhelixp[4],2))/m_lhelixp[2];
4362 m_lpt = 1/m_lhelix[2];
4363 m_lpte = 1/m_lhelixe[2];
4364 m_lptmu = 1/m_lhelixmu[2];
4365 m_lptk = 1/m_lhelixk[2];
4366 m_lptp = 1/m_lhelixp[2];
4368 m_fpt = 1/m_fhelix[2];
4369 m_fpte = 1/m_fhelixe[2];
4370 m_fptmu = 1/m_fhelixmu[2];
4371 m_fptk = 1/m_fhelixk[2];
4372 m_fptp = 1/m_fhelixp[2];
4375 std::cout<<
" "<<std::endl;
4376 std::cout<<
"in file Kalman_fitting_anal ,the m_fpt is .."<<m_fpt<<std::endl;
4377 std::cout<<
"in file Kalman_fitting_anal ,the m_fpte is .."<<m_fpte<<std::endl;
4378 std::cout<<
"in file Kalman_fitting_anal ,the m_fptmu is .."<<m_fptmu<<std::endl;
4379 std::cout<<
"in file Kalman_fitting_anal ,the m_fptk is .."<<m_fptk<<std::endl;
4380 std::cout<<
"in file Kalman_fitting_anal ,the m_fptp is .."<<m_fptp<<std::endl;
4383 m_zpt = 1/m_zhelix[2];
4384 m_zpte = 1/m_zhelixe[2];
4385 m_zptmu = 1/m_zhelixmu[2];
4386 m_zptk = 1/m_zhelixk[2];
4387 m_zptp = 1/m_zhelixp[2];
4390 std::cout<<
"in file Kalman_fitting_anal ,the m_zpt is .."<<m_zpt<<std::endl;
4391 std::cout<<
"in file Kalman_fitting_anal ,the m_zpte is .."<<m_zpte<<std::endl;
4392 std::cout<<
"in file Kalman_fitting_anal ,the m_zptmu is .."<<m_zptmu<<std::endl;
4393 std::cout<<
"in file Kalman_fitting_anal ,the m_zptk is .."<<m_zptk<<std::endl;
4394 std::cout<<
"in file Kalman_fitting_anal ,the m_zptp is .."<<m_zptp<<std::endl;
4396 m_zptot = sqrt(1+pow(m_zhelix[4],2))/m_zhelix[2];
4397 m_zptote = sqrt(1+pow(m_zhelixe[4],2))/m_zhelixe[2];
4398 m_zptotmu = sqrt(1+pow(m_zhelixmu[4],2))/m_zhelixmu[2];
4399 m_zptotk = sqrt(1+pow(m_zhelixk[4],2))/m_zhelixk[2];
4400 m_zptotp = sqrt(1+pow(m_zhelixp[4],2))/m_zhelixp[2];
4403 std::cout<<
"in file Kalman_fitting_anal ,the m_zptot is .."<<m_zptot<<std::endl;
4404 std::cout<<
"in file Kalman_fitting_anal ,the m_zptote is .."<<m_zptote<<std::endl;
4405 std::cout<<
"in file Kalman_fitting_anal ,the m_zptotmu is .."<<m_zptotmu<<std::endl;
4406 std::cout<<
"in file Kalman_fitting_anal ,the m_zptotk is .."<<m_zptotk<<std::endl;
4407 std::cout<<
"in file Kalman_fitting_anal ,the m_zptotp is .."<<m_zptotp<<std::endl;
4411 m_zsigp = sqrt(pow((m_zptot/m_zhelix[2]),2)*m_zerror[5]+
4412 pow((m_zhelix[4]/m_zptot),2)*pow((1/m_zhelix[2]),4)*m_zerror[14]-
4413 2*m_zhelix[4]*m_zerror[12]*pow((1/m_zhelix[2]),3));
4414 m_zsigpe = sqrt(pow((m_zptote/m_zhelixe[2]),2)*m_zerrore[5]+
4415 pow((m_zhelixe[4]/m_zptote),2)*pow((1/m_zhelixe[2]),4)*m_zerrore[14]-
4416 2*m_zhelixe[4]*m_zerrore[12]*pow((1/m_zhelixe[2]),3));
4417 m_zsigpmu = sqrt(pow((m_zptotmu/m_zhelixmu[2]),2)*m_zerrormu[5]+
4418 pow((m_zhelixmu[4]/m_zptotmu),2)*pow((1/m_zhelixmu[2]),4)*m_zerrormu[14]-
4419 2*m_zhelixmu[4]*m_zerrormu[12]*pow((1/m_zhelixmu[2]),3));
4420 m_zsigpk = sqrt(pow((m_zptotk/m_zhelixk[2]),2)*m_zerrork[5]+
4421 pow((m_zhelixk[4]/m_zptotk),2)*pow((1/m_zhelixk[2]),4)*m_zerrork[14]-
4422 2*m_zhelixk[4]*m_zerrork[12]*pow((1/m_zhelixk[2]),3));
4423 m_zsigpp = sqrt(pow((m_zptotp/m_zhelixp[2]),2)*m_zerrorp[5]+
4424 pow((m_zhelixp[4]/m_zptotp),2)*pow((1/m_zhelixp[2]),4)*m_zerrorp[14]-
4425 2*m_zhelixp[4]*m_zerrorp[12]*pow((1/m_zhelixp[2]),3));
4428 StatusCode sc1 = m_nt1->write();
4429 if( sc1.isFailure() ) cout<<
"Ntuple1 filling failed!"<<endl;
4434 phi1 = (*iter_trk)->getFFi0();
4435 r1 = (*iter_trk)->getFDr();
4436 z1 = (*iter_trk)->getFDz();
4437 kap1 = (*iter_trk)->getFCpa();
4438 tanl1 = (*iter_trk)->getFTanl();
4439 charge1 = kap1/fabs(kap1);
4442 p1 = sqrt(1+tanl1*tanl1)/kap1;
4443 the1 =
M_PI/2-atan(tanl1);
4446 pz1= tanl1/fabs(kap1);
4448 }
else if(jj == 2) {
4449 phi2 = (*iter_trk)->getFFi0();
4450 r2 = (*iter_trk)->getFDr();
4451 z2 = (*iter_trk)->getFDz();
4452 kap2 = (*iter_trk)->getFCpa();
4453 tanl2 = (*iter_trk)->getFTanl();
4454 charge2 = kap2/fabs(kap2);
4457 p2 = sqrt(1+tanl2*tanl2)/kap1;
4458 the2 =
M_PI/2-atan(tanl2);
4461 pz2= tanl2/fabs(kap2);
4469 m_delthe = the1 + the2;
4472 m_delpx = charge1*fabs(px1) + charge2*fabs(px2);
4473 m_delpy = charge1*fabs(py1) + charge2*fabs(py2);
4474 m_delpz = charge1*fabs(pz1) + charge2*fabs(pz2);
4476 StatusCode sc2 = m_nt2->write();
4477 if( sc2.isFailure() ) cout<<
"Ntuple2 filling failed!"<<endl;
4486 cout <<
"Kalfitting finished " << std::endl;
4493 MsgStream log(
msgSvc(), name());
4494 double Pt_threshold(0.3);
4495 Hep3Vector IP(0,0,0);
4502 if ( !&whMgr )
return;
4505 int ntrk = mdcMgr->size();
4506 double* rPt =
new double[ntrk];
4507 int* rOM =
new int[ntrk];
4508 unsigned int* order =
new unsigned int[ntrk];
4509 unsigned int* rCont =
new unsigned int[ntrk];
4510 unsigned int* rGen =
new unsigned int[ntrk];
4513 for(vector<MdcRec_trk>::iterator it = mdcMgr->begin(),
4514 end = mdcMgr->end(); it != end; it++) {
4518 rPt[index] = 1 / fabs(it->helix[2]);
4519 if(
debug_ == 4) cout<<
"rPt...[ "<<index<<
" ]...."<< rPt[index] <<endl;
4520 if(rPt[index] < 0) rPt[index] =
DBL_MAX;
4522 std::vector<MdcRec_wirhit*> pt = it->hitcol ;
4523 if(
debug_ == 4) cout<<
"ppt size: "<< pt.size()<<endl;
4525 for (vector<MdcRec_wirhit*>::iterator ii = pt.end()-1;
4526 ii !=pt.begin()-1; ii--) {
4527 int lyr((*ii)->geo->Lyr()->Id());
4528 if (outermost < lyr) outermost = lyr;
4529 if(
debug_ == 4) cout<<
"outmost: "<<outermost<<
" lyr: "<<lyr<<endl;
4531 rOM[index] = outermost;
4532 order[index] = index;
4537 for (
int j, k = ntrk - 1; k >= 0; k = j){
4539 for(
int i = 1; i <= k; i++)
4540 if(rPt[i - 1] < rPt[i]){
4542 std::swap(order[i], order[j]);
4543 std::swap(rPt[i], rPt[j]);
4544 std::swap(rOM[i], rOM[j]);
4545 std::swap(rCont[i], rCont[j]);
4546 std::swap(rGen[i], rGen[j]);
4553 DataObject *aReconEvent;
4554 eventSvc()->findObject(
"/Event/Recon",aReconEvent);
4558 StatusCode sc = eventSvc()->registerObject(
"/Event/Recon",recevt );
4559 if(sc!=StatusCode::SUCCESS) {
4560 log << MSG::FATAL <<
"Could not register ReconEvent" <<endreq;
4568 log << MSG::INFO <<
"beginning to make RecMdcKalTrackCol" <<endreq;
4571 for(
int l = 0; l < ntrk; l++) {
4573 MdcRec_trk& TrasanTRK = *(mdcMgr->begin() + order[l]);
4577 int trasqual = TrasanTRK_add.
quality;
4578 if(
debug_ == 4) cout<<
"kalman_fitting>trasqual... "<<trasqual<<endl;
4579 if (trasqual)
continue;
4583 cout <<
"******* KalFit NUMBER : " << newcount << std::endl;
4587 if ((TrasanTRK_add.
decision & 32) == 32 ||
4588 (TrasanTRK_add.
decision & 64) == 64) type = 1;
4593 TrasanTRK.
pivot[2]);
4596 for(
int i = 0; i < 5; i++)
4597 a[i] = TrasanTRK.
helix[i];
4600 for(
int i = 0, k = 0; i < 5; i++) {
4601 for(
int j = 0; j <= i; j++) {
4603 ea[j][i] = ea[i][j];
4609 double fiTerm = TrasanTRK.
fiTerm;
4615 int inlyr(999), outlyr(-1);
4616 int* rStat =
new int[43];
4617 for(
int irStat=0;irStat<43;++irStat)rStat[irStat]=0;
4618 std::vector<MdcRec_wirhit*> pt=TrasanTRK.
hitcol;
4620 if(
debug_ == 4) cout<<
"*********Pt size****"<< pt.size()<<endl;
4623 for (vector<MdcRec_wirhit*>::iterator ii = pt.end()-1;
4624 ii != pt.begin()-1; ii--) {
4625 Num[(*ii)->geo->Lyr()->Id()]++;
4629 for (vector<MdcRec_wirhit*>::iterator ii = pt.end()-1;
4630 ii != pt.begin()-1; ii--) {
4633 if (Num[(*ii)->geo->Lyr()->Id()]>3) {
4635 cout <<
"WARNING: I found " << Num[(*ii)->geo->Lyr()->Id()]
4636 <<
" hits in the layer "
4637 << (*ii)->geo->Lyr()->Id() << std::endl;
4664 double dist[2] = {rechit.
ddl, rechit.
ddr};
4665 double erdist[2] = {rechit.
erddl, rechit.
erddr};
4670 if (rechit.
lr==2 || rechit.
lr==0) lr_decision=-1;
4672 else if (rechit.
lr==1) lr_decision=1;
4675 int ind(geo->
Lyr()->
Id());
4677 lr_decision, rechit.
tdc,
4682 if (inlyr>ind) inlyr = ind;
4683 if (outlyr<ind) outlyr = ind;
4686 cout <<
"**** NUMBER OF Mdc HITS (TRASAN) = " << hit_asso << std::endl;
4689 int empty_between(0), empty(0);
4690 for (
int i= inlyr; i <= outlyr; i++)
4691 if (!rStat[i]) empty_between++;
4692 empty = empty_between+inlyr+(42-outlyr);
4697 track_lead.
type(type);
4698 unsigned int nhit = track_lead.
HitsMdc().size();
4699 if (!nhit &&
debug_ == 4) {
4700 cout <<
" ATTENTION TRACK WITH ONLY HITS " << nhit << std::endl;
4705 double KalFitst(0), KalFitax(0), KalFitschi2(0);
4707 Hep3Vector outer_pivot(track_lead.
x(fiTerm));
4710 std::cout<<
"before track_lead.pivot(outer_pivot) ,the error matrix of track_lead is .."<<track_lead.
Ea()<<std::endl;
4712 track_lead.
pivot(outer_pivot);
4717 HepSymMatrix Eakal(5,0);
4721 double costheta = track_lead.
a()[4] / sqrt(1.0 + track_lead.
a()[4]*track_lead.
a()[4]);
4731 std::cout <<
"from Mdc Pattern Recognition: " << std::endl;
4737 std::cout <<
" dr = " << work.
a()[0]
4738 <<
", Er_dr = " << sqrt(work.
Ea()[0][0]) << std::endl;
4739 std::cout <<
" phi0 = " << work.
a()[1]
4740 <<
", Er_phi0 = " << sqrt(work.
Ea()[1][1]) << std::endl;
4741 std::cout <<
" PT = " << 1/work.
a()[2]
4742 <<
", Er_kappa = " << sqrt(work.
Ea()[2][2]) << std::endl;
4743 std::cout <<
" dz = " << work.
a()[3]
4744 <<
", Er_dz = " << sqrt(work.
Ea()[3][3]) << std::endl;
4745 std::cout <<
" tanl = " << work.
a()[4]
4746 <<
", Er_tanl = " << sqrt(work.
Ea()[4][4]) << std::endl;
4754 cout <<
" Mdc FIRST KALMAN FIT " << std::endl;
4759 cout <<
" dr = " << work.
a()[0]
4760 <<
", Er_dr = " << sqrt(work.
Ea()[0][0]) << std::endl;
4761 cout <<
" phi0 = " << work.
a()[1]
4762 <<
", Er_phi0 = " << sqrt(work.
Ea()[1][1]) << std::endl;
4763 cout <<
" PT = " << 1/work.
a()[2]
4764 <<
", Er_kappa = " << sqrt(work.
Ea()[2][2]) << std::endl;
4765 cout <<
" dz = " << work.
a()[3]
4766 <<
", Er_dz = " << sqrt(work.
Ea()[3][3]) << std::endl;
4767 cout <<
" tanl = " << work.
a()[4]
4768 <<
", Er_tanl = " << sqrt(work.
Ea()[4][4]) << std::endl;
4775 complete_track(TrasanTRK, TrasanTRK_add, track_lead, kaltrk,kalcol,segcol);
4781 DataObject *aRecKalEvent;
4782 eventSvc()->findObject(
"/Event/Recon/RecMdcKalTrackCol", aRecKalEvent);
4783 if(aRecKalEvent!=NULL) {
4785 kalsc = eventSvc()->unregisterObject(
"/Event/Recon/RecMdcKalTrackCol");
4786 if(kalsc != StatusCode::SUCCESS) {
4787 log << MSG::FATAL <<
"Could not unregister RecMdcKalTrack collection" << endreq;
4792 kalsc = eventSvc()->registerObject(
"/Event/Recon/RecMdcKalTrackCol", kalcol);
4793 if( kalsc.isFailure()) {
4794 log << MSG::FATAL <<
"Could not register RecMdcKalTrack" << endreq;
4797 log << MSG::INFO <<
"RecMdcKalTrackCol registered successfully!" <<endreq;
4803 DataObject *aRecKalSegEvent;
4804 eventSvc()->findObject(
"/Event/Recon/RecMdcKalHelixSegCol", aRecKalSegEvent);
4805 if(aRecKalSegEvent!=NULL) {
4807 segsc = eventSvc()->unregisterObject(
"/Event/Recon/RecMdcKalHelixSegCol");
4808 if(segsc != StatusCode::SUCCESS) {
4809 log << MSG::FATAL <<
"Could not unregister RecMdcKalHelixSegCol collection" << endreq;
4814 segsc = eventSvc()->registerObject(
"/Event/Recon/RecMdcKalHelixSegCol", segcol);
4815 if( segsc.isFailure() ) {
4816 log << MSG::FATAL <<
"Could not register RecMdcKalHelixSeg" << endreq;
4819 log << MSG::INFO <<
"RecMdcKalHelixSegCol registered successfully!" <<endreq;
4822 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.);
4823 double r1(0.),r2(0.),kap1(999.),kap2(999.),tanl1(0.),tanl2(0.);
4826 SmartDataPtr<RecMdcKalTrackCol> kaltrkCol(eventSvc(),
"/Event/Recon/RecMdcKalTrackCol");
4828 log << MSG::FATAL <<
"Could not find RecMdcKalTrackCol" << endreq;
4831 log << MSG::INFO <<
"Begin to check RecMdcKalTrackCol"<<endreq;
4832 RecMdcKalTrackCol::iterator iter_trk = kaltrkCol->begin();
4833 for(
int jj=1; iter_trk != kaltrkCol->end(); iter_trk++,jj++) {
4834 log << MSG::DEBUG <<
"retrieved MDC Kalmantrack:"
4835 <<
"Track Id: " << (*iter_trk)->getTrackId()
4836 <<
" Mass of the fit: "<< (*iter_trk)->getMass(2)<< endreq
4837 <<
" Length of the track: "<< (*iter_trk)->getLength(2)
4838 <<
" Tof of the track: "<< (*iter_trk)->getTof(2) << endreq
4839 <<
" Chisq of the fit: "<< (*iter_trk)->getChisq(0,2)
4840 <<
" "<< (*iter_trk)->getChisq(1,2) << endreq
4841 <<
"Ndf of the fit: "<< (*iter_trk)->getNdf(0,1)
4842 <<
" "<< (*iter_trk)->getNdf(1,2) << endreq
4843 <<
"Kappa " << (*iter_trk)->getZHelix()[2]
4848 std::cout<<
"the size of gothelixsegs ..."<<gothelixsegs.size()<<std::endl;
4851 HelixSegRefVec::iterator it_gothelixseg = gothelixsegs.begin();
4852 for( ; it_gothelixseg != gothelixsegs.end(); it_gothelixseg++) {
4854 std::cout<<
"the layerId of this helixseg is ..."<<(*it_gothelixseg)->getLayerId()<<std::endl;
4855 std::cout<<
"the residual of this helixseg exclude the meas hit"<<(*it_gothelixseg)->getResExcl()<<std::endl;
4856 std::cout<<
"the residual of this helixseg include the meas hit"<<(*it_gothelixseg)->getResIncl()<<std::endl;
4857 std::cout<<
"the track id of the helixseg is ..."<<(*it_gothelixseg)->getTrackId() <<std::endl;
4858 std::cout<<
"the tof of the helixseg is ..."<<(*it_gothelixseg)->getTof()<<std::endl;
4859 std::cout<<
"the Zhit of the helixseg is ..."<<(*it_gothelixseg)->getZhit()<<std::endl;
4862 for(
int i = 0; i<43; i++) {
4863 log << MSG::DEBUG <<
"retrieved pathl["<<i<<
"]= "
4864 << (*iter_trk)->getPathl(i) <<endreq;
4868 m_trackid = (*iter_trk)->getTrackId();
4869 for(
int jj =0, iii=0; jj<5; jj++) {
4870 m_length[jj] = (*iter_trk)->getLength(jj);
4871 m_tof[jj] = (*iter_trk)->getTof(jj);
4872 m_nhits[jj] = (*iter_trk)->getNhits(jj);
4873 m_zhelix[jj] = (*iter_trk)->getZHelix()[jj];
4874 m_zhelixe[jj] = (*iter_trk)->getZHelixE()[jj];
4875 m_zhelixmu[jj] = (*iter_trk)->getZHelixMu()[jj];
4876 m_zhelixk[jj] = (*iter_trk)->getZHelixK()[jj];
4877 m_zhelixp[jj] = (*iter_trk)->getZHelixP()[jj];
4878 m_fhelix[jj] = (*iter_trk)->getFHelix()[jj];
4879 m_fhelixe[jj] = (*iter_trk)->getFHelixE()[jj];
4880 m_fhelixmu[jj] = (*iter_trk)->getFHelixMu()[jj];
4881 m_fhelixk[jj] = (*iter_trk)->getFHelixK()[jj];
4882 m_fhelixp[jj] = (*iter_trk)->getFHelixP()[jj];
4883 m_lhelix[jj] = (*iter_trk)->getLHelix()[jj];
4884 m_lhelixe[jj] = (*iter_trk)->getLHelixE()[jj];
4885 m_lhelixmu[jj] = (*iter_trk)->getLHelixMu()[jj];
4886 m_lhelixk[jj] = (*iter_trk)->getLHelixK()[jj];
4887 m_lhelixp[jj] = (*iter_trk)->getLHelixP()[jj];
4889 for(
int kk=0; kk<=jj; kk++,iii++) {
4890 m_zerror[iii] = (*iter_trk)->getZError()[jj][kk];
4891 m_zerrore[iii] = (*iter_trk)->getZErrorE()[jj][kk];
4892 m_zerrormu[iii] = (*iter_trk)->getZErrorMu()[jj][kk];
4893 m_zerrork[iii] = (*iter_trk)->getZErrorK()[jj][kk];
4894 m_zerrorp[iii] = (*iter_trk)->getZErrorP()[jj][kk];
4895 m_ferror[iii] = (*iter_trk)->getFError()[jj][kk];
4896 m_ferrore[iii] = (*iter_trk)->getFErrorE()[jj][kk];
4897 m_ferrormu[iii] = (*iter_trk)->getFErrorMu()[jj][kk];
4898 m_ferrork[iii] = (*iter_trk)->getFErrorK()[jj][kk];
4899 m_ferrorp[iii] = (*iter_trk)->getFErrorP()[jj][kk];
4900 m_lerror[iii] = (*iter_trk)->getLError()[jj][kk];
4901 m_lerrore[iii] = (*iter_trk)->getLErrorE()[jj][kk];
4902 m_lerrormu[iii] = (*iter_trk)->getLErrorMu()[jj][kk];
4903 m_lerrork[iii] = (*iter_trk)->getLErrorK()[jj][kk];
4904 m_lerrorp[iii] = (*iter_trk)->getLErrorP()[jj][kk];
4945 m_chisq[0][0] = (*iter_trk)->getChisq(0,0);
4946 m_chisq[1][0] = (*iter_trk)->getChisq(0,1);
4947 m_chisq[2][0] = (*iter_trk)->getChisq(0,2);
4948 m_chisq[3][0] = (*iter_trk)->getChisq(0,3);
4949 m_chisq[4][0] = (*iter_trk)->getChisq(0,4);
4950 m_chisq[0][1] = (*iter_trk)->getChisq(1,0);
4951 m_chisq[1][1] = (*iter_trk)->getChisq(1,1);
4952 m_chisq[2][1] = (*iter_trk)->getChisq(1,2);
4953 m_chisq[3][1] = (*iter_trk)->getChisq(1,3);
4954 m_chisq[4][1] = (*iter_trk)->getChisq(1,4);
4956 m_ndf[0][0] = (*iter_trk)->getNdf(0,0);
4957 m_ndf[1][0] = (*iter_trk)->getNdf(0,1);
4958 m_ndf[2][0] = (*iter_trk)->getNdf(0,2);
4959 m_ndf[3][0] = (*iter_trk)->getNdf(0,3);
4960 m_ndf[4][0] = (*iter_trk)->getNdf(0,4);
4961 m_ndf[0][1] = (*iter_trk)->getNdf(1,0);
4962 m_ndf[1][1] = (*iter_trk)->getNdf(1,1);
4963 m_ndf[2][1] = (*iter_trk)->getNdf(1,2);
4964 m_ndf[3][1] = (*iter_trk)->getNdf(1,3);
4965 m_ndf[4][1] = (*iter_trk)->getNdf(1,4);
4967 m_stat[0][0] = (*iter_trk)->getStat(0,0);
4968 m_stat[1][0] = (*iter_trk)->getStat(0,1);
4969 m_stat[2][0] = (*iter_trk)->getStat(0,2);
4970 m_stat[3][0] = (*iter_trk)->getStat(0,3);
4971 m_stat[4][0] = (*iter_trk)->getStat(0,4);
4972 m_stat[0][1] = (*iter_trk)->getStat(1,0);
4973 m_stat[1][1] = (*iter_trk)->getStat(1,1);
4974 m_stat[2][1] = (*iter_trk)->getStat(1,2);
4975 m_stat[3][1] = (*iter_trk)->getStat(1,3);
4976 m_stat[4][1] = (*iter_trk)->getStat(1,4);
4978 m_fptot = sqrt(1+pow(m_fhelix[4],2))/m_fhelix[2];
4979 m_fptote = sqrt(1+pow(m_fhelixe[4],2))/m_fhelixe[2];
4980 m_fptotmu = sqrt(1+pow(m_fhelixmu[4],2))/m_fhelixmu[2];
4981 m_fptotk = sqrt(1+pow(m_fhelixk[4],2))/m_fhelixk[2];
4982 m_fptotp = sqrt(1+pow(m_fhelixp[4],2))/m_fhelixp[2];
4984 m_zpt = 1/m_zhelix[2];
4985 m_zpte = 1/m_zhelixe[2];
4986 m_zptmu = 1/m_zhelixmu[2];
4987 m_zptk = 1/m_zhelixk[2];
4988 m_zptp = 1/m_zhelixp[2];
4990 m_fpt = 1/m_fhelix[2];
4991 m_fpte = 1/m_fhelixe[2];
4992 m_fptmu = 1/m_fhelixmu[2];
4993 m_fptk = 1/m_fhelixk[2];
4994 m_fptp = 1/m_fhelixp[2];
4996 m_lpt = 1/m_lhelix[2];
4997 m_lpte = 1/m_lhelixe[2];
4998 m_lptmu = 1/m_lhelixmu[2];
4999 m_lptk = 1/m_lhelixk[2];
5000 m_lptp = 1/m_lhelixp[2];
5002 m_lptot = sqrt(1+pow(m_lhelix[4],2))/m_lhelix[2];
5003 m_lptote = sqrt(1+pow(m_lhelixe[4],2))/m_lhelixe[2];
5004 m_lptotmu = sqrt(1+pow(m_lhelixmu[4],2))/m_lhelixmu[2];
5005 m_lptotk = sqrt(1+pow(m_lhelixk[4],2))/m_lhelixk[2];
5006 m_lptotp = sqrt(1+pow(m_lhelixp[4],2))/m_lhelixp[2];
5008 m_zptot = sqrt(1+pow(m_zhelix[4],2))/m_zhelix[2];
5009 m_zptote = sqrt(1+pow(m_zhelixe[4],2))/m_zhelixe[2];
5010 m_zptotmu = sqrt(1+pow(m_zhelixmu[4],2))/m_zhelixmu[2];
5011 m_zptotk = sqrt(1+pow(m_zhelixk[4],2))/m_zhelixk[2];
5012 m_zptotp = sqrt(1+pow(m_zhelixp[4],2))/m_zhelixp[2];
5014 m_zsigp = sqrt(pow((m_zptot/m_zhelix[2]),2)*m_zerror[5]+
5015 pow((m_zhelix[4]/m_zptot),2)*pow((1/m_zhelix[2]),4)*m_zerror[14]-
5016 2*m_zhelix[4]*m_zerror[12]*pow((1/m_zhelix[2]),3));
5017 m_zsigpe = sqrt(pow((m_zptote/m_zhelixe[2]),2)*m_zerrore[5]+
5018 pow((m_zhelixe[4]/m_zptote),2)*pow((1/m_zhelixe[2]),4)*m_zerrore[14]-
5019 2*m_zhelixe[4]*m_zerrore[12]*pow((1/m_zhelixe[2]),3));
5020 m_zsigpmu = sqrt(pow((m_zptotmu/m_zhelixmu[2]),2)*m_zerrormu[5]+
5021 pow((m_zhelixmu[4]/m_zptotmu),2)*pow((1/m_zhelixmu[2]),4)*m_zerrormu[14]-
5022 2*m_zhelixmu[4]*m_zerrormu[12]*pow((1/m_zhelixmu[2]),3));
5023 m_zsigpk = sqrt(pow((m_zptotk/m_zhelixk[2]),2)*m_zerrork[5]+
5024 pow((m_zhelixk[4]/m_zptotk),2)*pow((1/m_zhelixk[2]),4)*m_zerrork[14]-
5025 2*m_zhelixk[4]*m_zerrork[12]*pow((1/m_zhelixk[2]),3));
5026 m_zsigpp = sqrt(pow((m_zptotp/m_zhelixp[2]),2)*m_zerrorp[5]+
5027 pow((m_zhelixp[4]/m_zptotp),2)*pow((1/m_zhelixp[2]),4)*m_zerrorp[14]-
5028 2*m_zhelixp[4]*m_zerrorp[12]*pow((1/m_zhelixp[2]),3));
5031 StatusCode sc1 = m_nt1->write();
5032 if( sc1.isFailure() ) cout<<
"Ntuple1 filling failed!"<<endl;
5037 phi1 = (*iter_trk)->getFFi0();
5038 r1 = (*iter_trk)->getFDr();
5039 z1 = (*iter_trk)->getFDz();
5040 kap1 = (*iter_trk)->getFCpa();
5041 tanl1 = (*iter_trk)->getFTanl();
5044 p1 = sqrt(1+tanl1*tanl1)/kap1;
5045 the1 =
M_PI/2-atan(tanl1);
5046 }
else if(jj == 2) {
5047 phi2 = (*iter_trk)->getFFi0();
5048 r2 = (*iter_trk)->getFDr();
5049 z2 = (*iter_trk)->getFDz();
5050 kap2 = (*iter_trk)->getFCpa();
5051 tanl2 = (*iter_trk)->getFTanl();
5054 p2 = sqrt(1+tanl2*tanl2)/kap1;
5055 the2 =
M_PI/2-atan(tanl2);
5063 m_delthe = the1 + the2;
5066 StatusCode sc2 = m_nt2->write();
5067 if( sc2.isFailure() ) cout<<
"Ntuple2 filling failed!"<<endl;
5075 cout <<
"Kalfitting finished " << std::endl;
5082 MsgStream log(
msgSvc(), name());
5083 double Pt_threshold(0.3);
5084 Hep3Vector IP(0,0,0);
5091 if ( !&whMgr )
return;
5094 int ntrk = mdcMgr->size();
5097 int nhits = whMgr->size();
5102 DataObject *aReconEvent;
5103 eventSvc()->findObject(
"/Event/Recon",aReconEvent);
5107 StatusCode sc = eventSvc()->registerObject(
"/Event/Recon",recevt );
5108 if(sc!=StatusCode::SUCCESS) {
5109 log << MSG::FATAL <<
"Could not register ReconEvent" <<endreq;
5117 log << MSG::INFO <<
"beginning to make RecMdcKalTrackCol" <<endreq;
5129 if ((TrasanTRK_add.
decision & 32) == 32 ||
5130 (TrasanTRK_add.
decision & 64) == 64) type = 1;
5135 TrasanTRK.
pivot[2]);
5138 for(
int i = 0; i < 5; i++)
5139 a[i] = TrasanTRK.
helix[i];
5142 for(
int i = 0, k = 0; i < 5; i++) {
5143 for(
int j = 0; j <= i; j++) {
5145 ea[j][i] = ea[i][j];
5151 double fiTerm = TrasanTRK.
fiTerm;
5159 int trasqual = TrasanTRK_add.
quality;
5160 if(
debug_ == 4) cout<<
"kalman_fitting>trasqual... "<<trasqual<<endl;
5161 if (trasqual)
return;
5163 int inlyr(999), outlyr(-1);
5164 int* rStat =
new int[43];
5165 for(
int irStat=0;irStat<43;++irStat)rStat[irStat]=0;
5166 std::vector<MdcRec_wirhit*> pt=TrasanTRK.
hitcol;
5168 if(
debug_ == 4) cout<<
"*********Pt size****"<< pt.size()<<endl;
5171 for (vector<MdcRec_wirhit*>::iterator ii = pt.end()-1;
5172 ii != pt.begin()-1; ii--) {
5173 Num[(*ii)->geo->Lyr()->Id()]++;
5176 for (vector<MdcRec_wirhit*>::iterator ii = pt.end()-1;
5177 ii != pt.begin()-1; ii--) {
5180 if (Num[(*ii)->geo->Lyr()->Id()]>3) {
5182 cout <<
"WARNING: I found " << Num[(*ii)->geo->Lyr()->Id()]
5183 <<
" hits in the layer "
5184 << (*ii)->geo->Lyr()->Id() << std::endl;
5190 double dist[2] = {rechit.
ddl, rechit.
ddr};
5191 double erdist[2] = {rechit.
erddl, rechit.
erddr};
5196 if (rechit.
lr==2 || rechit.
lr==0) lr_decision=-1;
5198 else if (rechit.
lr==1) lr_decision=1;
5201 int ind(geo->
Lyr()->
Id());
5203 lr_decision, rechit.
tdc,
5208 if (inlyr>ind) inlyr = ind;
5209 if (outlyr<ind) outlyr = ind;
5212 int empty_between(0), empty(0);
5213 for (
int i= inlyr; i <= outlyr; i++)
5214 if (!rStat[i]) empty_between++;
5215 empty = empty_between+inlyr+(42-outlyr);
5219 cout <<
"**** NUMBER OF Mdc HITS (TRASAN) = " << hit_asso << std::endl;
5224 track_lead.
type(type);
5225 unsigned int nhit = track_lead.
HitsMdc().size();
5227 cout <<
" ATTENTION TRACK WITH ONLY HITS " << nhit << std::endl;
5232 double KalFitst(0), KalFitax(0), KalFitschi2(0);
5234 Hep3Vector outer_pivot(track_lead.
x(fiTerm));
5237 std::cout<<
"before track_lead.pivot(outer_pivot) ,the error matrix of track_lead is .."<<track_lead.
Ea()<<std::endl;
5239 track_lead.
pivot(outer_pivot);
5244 HepSymMatrix Eakal(5,0);
5248 double costheta = track_lead.
a()[4] / sqrt(1.0 + track_lead.
a()[4]*track_lead.
a()[4]);
5258 std::cout <<
"from Mdc Pattern Recognition: " << std::endl;
5264 std::cout <<
" dr = " << work.
a()[0]
5265 <<
", Er_dr = " << sqrt(work.
Ea()[0][0]) << std::endl;
5266 std::cout <<
" phi0 = " << work.
a()[1]
5267 <<
", Er_phi0 = " << sqrt(work.
Ea()[1][1]) << std::endl;
5268 std::cout <<
" PT = " << 1/work.
a()[2]
5269 <<
", Er_kappa = " << sqrt(work.
Ea()[2][2]) << std::endl;
5270 std::cout <<
" dz = " << work.
a()[3]
5271 <<
", Er_dz = " << sqrt(work.
Ea()[3][3]) << std::endl;
5272 std::cout <<
" tanl = " << work.
a()[4]
5273 <<
", Er_tanl = " << sqrt(work.
Ea()[4][4]) << std::endl;
5280 cout <<
" Mdc FIRST KALMAN FIT " << std::endl;
5285 cout <<
" dr = " << work1.
a()[0]
5286 <<
", Er_dr = " << sqrt(work1.
Ea()[0][0]) << std::endl;
5287 cout <<
" phi0 = " << work1.
a()[1]
5288 <<
", Er_phi0 = " << sqrt(work1.
Ea()[1][1]) << std::endl;
5289 cout <<
" PT = " << 1/work1.
a()[2]
5290 <<
", Er_kappa = " << sqrt(work1.
Ea()[2][2]) << std::endl;
5291 cout <<
" dz = " << work1.
a()[3]
5292 <<
", Er_dz = " << sqrt(work1.
Ea()[3][3]) << std::endl;
5293 cout <<
" tanl = " << work1.
a()[4]
5294 <<
", Er_tanl = " << sqrt(work1.
Ea()[4][4]) << std::endl;
5301 complete_track(TrasanTRK, TrasanTRK_add, track_lead, kaltrk,kalcol,segcol);
5306 DataObject *aRecKalEvent;
5307 eventSvc()->findObject(
"/Event/Recon/RecMdcKalTrackCol", aRecKalEvent);
5308 if(aRecKalEvent!=NULL) {
5310 kalsc = eventSvc()->unregisterObject(
"/Event/Recon/RecMdcKalTrackCol");
5311 if(kalsc != StatusCode::SUCCESS) {
5312 log << MSG::FATAL <<
"Could not unregister RecMdcKalTrack collection" << endreq;
5317 kalsc = eventSvc()->registerObject(
"/Event/Recon/RecMdcKalTrackCol", kalcol);
5318 if( kalsc.isFailure()) {
5319 log << MSG::FATAL <<
"Could not register RecMdcKalTrack" << endreq;
5322 log << MSG::INFO <<
"RecMdcKalTrackCol registered successfully!" <<endreq;
5328 DataObject *aRecKalSegEvent;
5329 eventSvc()->findObject(
"/Event/Recon/RecMdcKalHelixSegCol", aRecKalSegEvent);
5330 if(aRecKalSegEvent!=NULL) {
5332 segsc = eventSvc()->unregisterObject(
"/Event/Recon/RecMdcKalHelixSegCol");
5333 if(segsc != StatusCode::SUCCESS) {
5334 log << MSG::FATAL <<
"Could not unregister RecMdcKalHelixSegCol collection" << endreq;
5339 segsc = eventSvc()->registerObject(
"/Event/Recon/RecMdcKalHelixSegCol", segcol);
5340 if( segsc.isFailure() ) {
5341 log << MSG::FATAL <<
"Could not register RecMdcKalHelixSeg" << endreq;
5344 log << MSG::INFO <<
"RecMdcKalHelixSegCol registered successfully!" <<endreq;
5347 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.);
5348 double r1(0.),r2(0.),kap1(999.),kap2(999.),tanl1(0.),tanl2(0.);
5351 SmartDataPtr<RecMdcKalTrackCol> kaltrkCol(eventSvc(),
"/Event/Recon/RecMdcKalTrackCol");
5353 log << MSG::FATAL <<
"Could not find RecMdcKalTrackCol" << endreq;
5356 log << MSG::INFO <<
"Begin to check RecMdcKalTrackCol"<<endreq;
5357 RecMdcKalTrackCol::iterator iter_trk = kaltrkCol->begin();
5358 for(
int jj=1; iter_trk != kaltrkCol->end(); iter_trk++,jj++) {
5359 log << MSG::DEBUG <<
"retrieved MDC Kalmantrack:"
5360 <<
"Track Id: " << (*iter_trk)->getTrackId()
5361 <<
" Mass of the fit: "<< (*iter_trk)->getMass(2)<< endreq
5362 <<
" Length of the track: "<< (*iter_trk)->getLength(2)
5363 <<
" Tof of the track: "<< (*iter_trk)->getTof(2) << endreq
5364 <<
" Chisq of the fit: "<< (*iter_trk)->getChisq(0,2)
5365 <<
" "<< (*iter_trk)->getChisq(1,2) << endreq
5366 <<
"Ndf of the fit: "<< (*iter_trk)->getNdf(0,1)
5367 <<
" "<< (*iter_trk)->getNdf(1,2) << endreq
5368 <<
"Kappa " << (*iter_trk)->getZHelix()[2]
5369 <<
"zhelixmu "<<(*iter_trk)->getZHelixMu()
5374 std::cout<<
"the size of gothelixsegs ..."<<gothelixsegs.size()<<std::endl;
5377 HelixSegRefVec::iterator it_gothelixseg = gothelixsegs.begin();
5378 for( ; it_gothelixseg != gothelixsegs.end(); it_gothelixseg++) {
5380 std::cout<<
"the layerId of this helixseg is ..."<<(*it_gothelixseg)->getLayerId()<<std::endl;
5381 std::cout<<
"the residual of this helixseg exclude the meas hit"<<(*it_gothelixseg)->getResExcl()<<std::endl;
5382 std::cout<<
"the residual of this helixseg include the meas hit"<<(*it_gothelixseg)->getResIncl()<<std::endl;
5383 std::cout<<
"the track id of the helixseg is ..."<<(*it_gothelixseg)->getTrackId() <<std::endl;
5384 std::cout<<
"the tof of the helixseg is ..."<<(*it_gothelixseg)->getTof()<<std::endl;
5385 std::cout<<
"the Zhit of the helixseg is ..."<<(*it_gothelixseg)->getZhit()<<std::endl;
5388 for(
int i = 0; i<43; i++) {
5389 log << MSG::DEBUG <<
"retrieved pathl["<<i<<
"]= "
5390 << (*iter_trk)->getPathl(i) <<endreq;
5394 m_trackid = (*iter_trk)->getTrackId();
5395 for(
int jj =0, iii=0; jj<5; jj++) {
5396 m_length[jj] = (*iter_trk)->getLength(jj);
5397 m_tof[jj] = (*iter_trk)->getTof(jj);
5398 m_nhits[jj] = (*iter_trk)->getNhits(jj);
5399 m_zhelix[jj] = (*iter_trk)->getZHelix()[jj];
5400 m_zhelixe[jj] = (*iter_trk)->getZHelixE()[jj];
5401 m_zhelixmu[jj] = (*iter_trk)->getZHelixMu()[jj];
5402 m_zhelixk[jj] = (*iter_trk)->getZHelixK()[jj];
5403 m_zhelixp[jj] = (*iter_trk)->getZHelixP()[jj];
5404 m_fhelix[jj] = (*iter_trk)->getFHelix()[jj];
5405 m_fhelixe[jj] = (*iter_trk)->getFHelixE()[jj];
5406 m_fhelixmu[jj] = (*iter_trk)->getFHelixMu()[jj];
5407 m_fhelixk[jj] = (*iter_trk)->getFHelixK()[jj];
5408 m_fhelixp[jj] = (*iter_trk)->getFHelixP()[jj];
5409 m_lhelix[jj] = (*iter_trk)->getLHelix()[jj];
5410 m_lhelixe[jj] = (*iter_trk)->getLHelixE()[jj];
5411 m_lhelixmu[jj] = (*iter_trk)->getLHelixMu()[jj];
5412 m_lhelixk[jj] = (*iter_trk)->getLHelixK()[jj];
5413 m_lhelixp[jj] = (*iter_trk)->getLHelixP()[jj];
5415 for(
int kk=0; kk<=jj; kk++,iii++) {
5416 m_zerror[iii] = (*iter_trk)->getZError()[jj][kk];
5417 m_zerrore[iii] = (*iter_trk)->getZErrorE()[jj][kk];
5418 m_zerrormu[iii] = (*iter_trk)->getZErrorMu()[jj][kk];
5419 m_zerrork[iii] = (*iter_trk)->getZErrorK()[jj][kk];
5420 m_zerrorp[iii] = (*iter_trk)->getZErrorP()[jj][kk];
5421 m_ferror[iii] = (*iter_trk)->getFError()[jj][kk];
5422 m_ferrore[iii] = (*iter_trk)->getFErrorE()[jj][kk];
5423 m_ferrormu[iii] = (*iter_trk)->getFErrorMu()[jj][kk];
5424 m_ferrork[iii] = (*iter_trk)->getFErrorK()[jj][kk];
5425 m_ferrorp[iii] = (*iter_trk)->getFErrorP()[jj][kk];
5426 m_lerror[iii] = (*iter_trk)->getLError()[jj][kk];
5427 m_lerrore[iii] = (*iter_trk)->getLErrorE()[jj][kk];
5428 m_lerrormu[iii] = (*iter_trk)->getLErrorMu()[jj][kk];
5429 m_lerrork[iii] = (*iter_trk)->getLErrorK()[jj][kk];
5430 m_lerrorp[iii] = (*iter_trk)->getLErrorP()[jj][kk];
5437 m_chisq[0][0] = (*iter_trk)->getChisq(0,0);
5438 m_chisq[1][0] = (*iter_trk)->getChisq(0,1);
5439 m_chisq[2][0] = (*iter_trk)->getChisq(0,2);
5440 m_chisq[3][0] = (*iter_trk)->getChisq(0,3);
5441 m_chisq[4][0] = (*iter_trk)->getChisq(0,4);
5442 m_chisq[0][1] = (*iter_trk)->getChisq(1,0);
5443 m_chisq[1][1] = (*iter_trk)->getChisq(1,1);
5444 m_chisq[2][1] = (*iter_trk)->getChisq(1,2);
5445 m_chisq[3][1] = (*iter_trk)->getChisq(1,3);
5446 m_chisq[4][1] = (*iter_trk)->getChisq(1,4);
5448 m_ndf[0][0] = (*iter_trk)->getNdf(0,0);
5449 m_ndf[1][0] = (*iter_trk)->getNdf(0,1);
5450 m_ndf[2][0] = (*iter_trk)->getNdf(0,2);
5451 m_ndf[3][0] = (*iter_trk)->getNdf(0,3);
5452 m_ndf[4][0] = (*iter_trk)->getNdf(0,4);
5453 m_ndf[0][1] = (*iter_trk)->getNdf(1,0);
5454 m_ndf[1][1] = (*iter_trk)->getNdf(1,1);
5455 m_ndf[2][1] = (*iter_trk)->getNdf(1,2);
5456 m_ndf[3][1] = (*iter_trk)->getNdf(1,3);
5457 m_ndf[4][1] = (*iter_trk)->getNdf(1,4);
5459 m_stat[0][0] = (*iter_trk)->getStat(0,0);
5460 m_stat[1][0] = (*iter_trk)->getStat(0,1);
5461 m_stat[2][0] = (*iter_trk)->getStat(0,2);
5462 m_stat[3][0] = (*iter_trk)->getStat(0,3);
5463 m_stat[4][0] = (*iter_trk)->getStat(0,4);
5464 m_stat[0][1] = (*iter_trk)->getStat(1,0);
5465 m_stat[1][1] = (*iter_trk)->getStat(1,1);
5466 m_stat[2][1] = (*iter_trk)->getStat(1,2);
5467 m_stat[3][1] = (*iter_trk)->getStat(1,3);
5468 m_stat[4][1] = (*iter_trk)->getStat(1,4);
5470 m_fptot = sqrt(1+pow(m_fhelix[4],2))/m_fhelix[2];
5471 m_fptote = sqrt(1+pow(m_fhelixe[4],2))/m_fhelixe[2];
5472 m_fptotmu = sqrt(1+pow(m_fhelixmu[4],2))/m_fhelixmu[2];
5473 m_fptotk = sqrt(1+pow(m_fhelixk[4],2))/m_fhelixk[2];
5474 m_fptotp = sqrt(1+pow(m_fhelixp[4],2))/m_fhelixp[2];
5476 m_zpt = 1/m_zhelix[2];
5477 m_zpte = 1/m_zhelixe[2];
5478 m_zptmu = 1/m_zhelixmu[2];
5479 m_zptk = 1/m_zhelixk[2];
5480 m_zptp = 1/m_zhelixp[2];
5482 m_fpt = 1/m_fhelix[2];
5483 m_fpte = 1/m_fhelixe[2];
5484 m_fptmu = 1/m_fhelixmu[2];
5485 m_fptk = 1/m_fhelixk[2];
5486 m_fptp = 1/m_fhelixp[2];
5488 m_lpt = 1/m_lhelix[2];
5489 m_lpte = 1/m_lhelixe[2];
5490 m_lptmu = 1/m_lhelixmu[2];
5491 m_lptk = 1/m_lhelixk[2];
5492 m_lptp = 1/m_lhelixp[2];
5494 m_lptot = sqrt(1+pow(m_lhelix[4],2))/m_lhelix[2];
5495 m_lptote = sqrt(1+pow(m_lhelixe[4],2))/m_lhelixe[2];
5496 m_lptotmu = sqrt(1+pow(m_lhelixmu[4],2))/m_lhelixmu[2];
5497 m_lptotk = sqrt(1+pow(m_lhelixk[4],2))/m_lhelixk[2];
5498 m_lptotp = sqrt(1+pow(m_lhelixp[4],2))/m_lhelixp[2];
5500 m_zptot = sqrt(1+pow(m_zhelix[4],2))/m_zhelix[2];
5501 m_zptote = sqrt(1+pow(m_zhelixe[4],2))/m_zhelixe[2];
5502 m_zptotmu = sqrt(1+pow(m_zhelixmu[4],2))/m_zhelixmu[2];
5503 m_zptotk = sqrt(1+pow(m_zhelixk[4],2))/m_zhelixk[2];
5504 m_zptotp = sqrt(1+pow(m_zhelixp[4],2))/m_zhelixp[2];
5506 m_zsigp = sqrt(pow((m_zptot/m_zhelix[2]),2)*m_zerror[5]+
5507 pow((m_zhelix[4]/m_zptot),2)*pow((1/m_zhelix[2]),4)*m_zerror[14]-
5508 2*m_zhelix[4]*m_zerror[12]*pow((1/m_zhelix[2]),3));
5509 m_zsigpe = sqrt(pow((m_zptote/m_zhelixe[2]),2)*m_zerrore[5]+
5510 pow((m_zhelixe[4]/m_zptote),2)*pow((1/m_zhelixe[2]),4)*m_zerrore[14]-
5511 2*m_zhelixe[4]*m_zerrore[12]*pow((1/m_zhelixe[2]),3));
5512 m_zsigpmu = sqrt(pow((m_zptotmu/m_zhelixmu[2]),2)*m_zerrormu[5]+
5513 pow((m_zhelixmu[4]/m_zptotmu),2)*pow((1/m_zhelixmu[2]),4)*m_zerrormu[14]-
5514 2*m_zhelixmu[4]*m_zerrormu[12]*pow((1/m_zhelixmu[2]),3));
5515 m_zsigpk = sqrt(pow((m_zptotk/m_zhelixk[2]),2)*m_zerrork[5]+
5516 pow((m_zhelixk[4]/m_zptotk),2)*pow((1/m_zhelixk[2]),4)*m_zerrork[14]-
5517 2*m_zhelixk[4]*m_zerrork[12]*pow((1/m_zhelixk[2]),3));
5518 m_zsigpp = sqrt(pow((m_zptotp/m_zhelixp[2]),2)*m_zerrorp[5]+
5519 pow((m_zhelixp[4]/m_zptotp),2)*pow((1/m_zhelixp[2]),4)*m_zerrorp[14]-
5520 2*m_zhelixp[4]*m_zerrorp[12]*pow((1/m_zhelixp[2]),3));
5523 StatusCode sc1 = m_nt1->write();
5524 if( sc1.isFailure() ) cout<<
"Ntuple1 filling failed!"<<endl;
5529 phi1 = (*iter_trk)->getFFi0();
5530 r1 = (*iter_trk)->getFDr();
5531 z1 = (*iter_trk)->getFDz();
5532 kap1 = (*iter_trk)->getFCpa();
5533 tanl1 = (*iter_trk)->getFTanl();
5536 p1 = sqrt(1+tanl1*tanl1)/kap1;
5537 the1 =
M_PI/2-atan(tanl1);
5538 }
else if(jj == 2) {
5539 phi2 = (*iter_trk)->getFFi0();
5540 r2 = (*iter_trk)->getFDr();
5541 z2 = (*iter_trk)->getFDz();
5542 kap2 = (*iter_trk)->getFCpa();
5543 tanl2 = (*iter_trk)->getFTanl();
5546 p2 = sqrt(1+tanl2*tanl2)/kap1;
5547 the2 =
M_PI/2-atan(tanl2);
5555 m_delthe = the1 + the2;
5558 StatusCode sc2 = m_nt2->write();
5559 if( sc2.isFailure() ) cout<<
"Ntuple2 filling failed!"<<endl;
5563 cout <<
"Kalfitting finished " << std::endl;
5571 MsgStream log(
msgSvc(), name());
5572 double Pt_threshold(0.3);
5573 Hep3Vector IP(0,0,0);
5580 if ( !&whMgr )
return;
5583 int ntrk = mdcMgr->size();
5586 int nhits = whMgr->size();
5590 double* rY =
new double[ntrk];
5591 double* rfiTerm =
new double[ntrk];
5592 double* rPt =
new double[ntrk];
5593 int* rOM =
new int[ntrk];
5594 unsigned int* order =
new unsigned int[ntrk];
5595 unsigned int* rCont =
new unsigned int[ntrk];
5596 unsigned int* rGen =
new unsigned int[ntrk];
5599 Hep3Vector csmp3[2];
5600 for(vector<MdcRec_trk>::iterator it = mdcMgr->begin(),
5601 end = mdcMgr->end(); it != end; it++) {
5603 rfiTerm[index]=it->fiTerm;
5608 rPt[index] = 1 / fabs(it->helix[2]);
5609 if(
debug_ == 4) cout<<
"rPt...[ "<<index<<
" ]...."<< rPt[index] <<endl;
5610 if(rPt[index] < 0) rPt[index] =
DBL_MAX;
5612 std::vector<MdcRec_wirhit*> pt = it->hitcol ;
5613 if(
debug_ == 4) cout<<
"ppt size: "<< pt.size()<<endl;
5615 for (vector<MdcRec_wirhit*>::iterator ii = pt.end()-1;
5616 ii !=pt.begin()-1; ii--) {
5617 int lyr((*ii)->geo->Lyr()->Id());
5618 if (outermost < lyr) {
5620 rY[index] = (*ii)->geo->Forward().y();
5622 if(
debug_ == 4) cout<<
"outmost: "<<outermost<<
" lyr: "<<lyr<<endl;
5624 rOM[index] = outermost;
5625 order[index] = index;
5630 for (
int j, k = ntrk - 1; k >= 0; k = j){
5632 for(
int i = 1; i <= k; i++)
5633 if(rY[i - 1] < rY[i]){
5635 std::swap(order[i], order[j]);
5636 std::swap(rY[i], rY[j]);
5637 std::swap(rOM[i], rOM[j]);
5638 std::swap(rCont[i], rCont[j]);
5639 std::swap(rGen[i], rGen[j]);
5648 DataObject *aReconEvent;
5649 eventSvc()->findObject(
"/Event/Recon",aReconEvent);
5653 StatusCode sc = eventSvc()->registerObject(
"/Event/Recon",recevt );
5654 if(sc!=StatusCode::SUCCESS) {
5655 log << MSG::FATAL <<
"Could not register ReconEvent" <<endreq;
5663 log << MSG::INFO <<
"beginning to make RecMdcKalTrackCol" <<endreq;
5670 MdcRec_trk& TrasanTRK = *(mdcMgr->begin() + order[1]);
5679 cout <<
"******* KalFit NUMBER : " << newcount << std::endl;
5683 if ((TrasanTRK_add.
decision & 32) == 32 ||
5684 (TrasanTRK_add.
decision & 64) == 64) type = 1;
5689 TrasanTRK.
pivot[2]);
5692 for(
int i = 0; i < 5; i++)
5693 a[i] = TrasanTRK.
helix[i];
5696 for(
int i = 0, k = 0; i < 5; i++) {
5697 for(
int j = 0; j <= i; j++) {
5699 ea[j][i] = ea[i][j];
5705 double fiTerm = TrasanTRK.
fiTerm;
5712 for(
int l = 0; l < ntrk; l++) {
5713 MdcRec_trk& TrasanTRK1 = *(mdcMgr->begin() + order[l]);
5714 MdcRec_trk_add& TrasanTRK_add1 = *(mdc_addMgr->begin()+order[l]);
5716 int trasqual = TrasanTRK_add1.
quality;
5717 if(
debug_ == 4) cout<<
"kalman_fitting>trasqual... "<<trasqual<<endl;
5718 if (trasqual)
continue;
5720 int inlyr(999), outlyr(-1);
5721 int* rStat =
new int[43];
5722 for(
int irStat=0;irStat<43;++irStat)rStat[irStat]=0;
5723 std::vector<MdcRec_wirhit*> pt=TrasanTRK1.
hitcol;
5725 if(
debug_ == 4) cout<<
"*********Pt size****"<< pt.size()<<endl;
5728 for (vector<MdcRec_wirhit*>::iterator ii = pt.end()-1;
5729 ii != pt.begin()-1; ii--) {
5730 Num[(*ii)->geo->Lyr()->Id()]++;
5733 for (vector<MdcRec_wirhit*>::iterator ii = pt.end()-1;
5734 ii != pt.begin()-1; ii--) {
5737 if (Num[(*ii)->geo->Lyr()->Id()]>3) {
5739 cout <<
"WARNING: I found " << Num[(*ii)->geo->Lyr()->Id()]
5740 <<
" hits in the layer "
5741 << (*ii)->geo->Lyr()->Id() << std::endl;
5747 double dist[2] = {rechit.
ddl, rechit.
ddr};
5748 double erdist[2] = {rechit.
erddl, rechit.
erddr};
5753 if (rechit.
lr==2 || rechit.
lr==0) lr_decision=-1;
5755 else if (rechit.
lr==1) lr_decision=1;
5758 int ind(geo->
Lyr()->
Id());
5760 lr_decision, rechit.
tdc,
5765 if (inlyr>ind) inlyr = ind;
5766 if (outlyr<ind) outlyr = ind;
5769 int empty_between(0), empty(0);
5770 for (
int i= inlyr; i <= outlyr; i++)
5771 if (!rStat[i]) empty_between++;
5772 empty = empty_between+inlyr+(42-outlyr);
5776 cout <<
"**** NUMBER OF Mdc HITS (TRASAN) = " << hit_asso << std::endl;
5781 track_lead.
type(type);
5782 unsigned int nhit = track_lead.
HitsMdc().size();
5784 cout <<
" ATTENTION TRACK WITH ONLY HITS " << nhit << std::endl;
5789 double KalFitst(0), KalFitax(0), KalFitschi2(0);
5791 Hep3Vector outer_pivot(track_lead.
x(fiTerm));
5794 std::cout<<
"before track_lead.pivot(outer_pivot) ,the error matrix of track_lead is .."<<track_lead.
Ea()<<std::endl;
5796 track_lead.
pivot(outer_pivot);
5801 HepSymMatrix Eakal(5,0);
5805 double costheta = track_lead.
a()[4] / sqrt(1.0 + track_lead.
a()[4]*track_lead.
a()[4]);
5815 std::cout <<
"from Mdc Pattern Recognition: " << std::endl;
5821 std::cout <<
" dr = " << work.
a()[0]
5822 <<
", Er_dr = " << sqrt(work.
Ea()[0][0]) << std::endl;
5823 std::cout <<
" phi0 = " << work.
a()[1]
5824 <<
", Er_phi0 = " << sqrt(work.
Ea()[1][1]) << std::endl;
5825 std::cout <<
" PT = " << 1/work.
a()[2]
5826 <<
", Er_kappa = " << sqrt(work.
Ea()[2][2]) << std::endl;
5827 std::cout <<
" dz = " << work.
a()[3]
5828 <<
", Er_dz = " << sqrt(work.
Ea()[3][3]) << std::endl;
5829 std::cout <<
" tanl = " << work.
a()[4]
5830 <<
", Er_tanl = " << sqrt(work.
Ea()[4][4]) << std::endl;
5837 cout <<
" Mdc FIRST KALMAN FIT " << std::endl;
5842 cout <<
" dr = " << work1.
a()[0]
5843 <<
", Er_dr = " << sqrt(work1.
Ea()[0][0]) << std::endl;
5844 cout <<
" phi0 = " << work1.
a()[1]
5845 <<
", Er_phi0 = " << sqrt(work1.
Ea()[1][1]) << std::endl;
5846 cout <<
" PT = " << 1/work1.
a()[2]
5847 <<
", Er_kappa = " << sqrt(work1.
Ea()[2][2]) << std::endl;
5848 cout <<
" dz = " << work1.
a()[3]
5849 <<
", Er_dz = " << sqrt(work1.
Ea()[3][3]) << std::endl;
5850 cout <<
" tanl = " << work1.
a()[4]
5851 <<
", Er_tanl = " << sqrt(work1.
Ea()[4][4]) << std::endl;
5858 complete_track(TrasanTRK, TrasanTRK_add, track_lead, kaltrk,kalcol,segcol);
5864 DataObject *aRecKalEvent;
5865 eventSvc()->findObject(
"/Event/Recon/RecMdcKalTrackCol", aRecKalEvent);
5866 if(aRecKalEvent!=NULL) {
5868 kalsc = eventSvc()->unregisterObject(
"/Event/Recon/RecMdcKalTrackCol");
5869 if(kalsc != StatusCode::SUCCESS) {
5870 log << MSG::FATAL <<
"Could not unregister RecMdcKalTrack collection" << endreq;
5875 kalsc = eventSvc()->registerObject(
"/Event/Recon/RecMdcKalTrackCol", kalcol);
5876 if( kalsc.isFailure()) {
5877 log << MSG::FATAL <<
"Could not register RecMdcKalTrack" << endreq;
5880 log << MSG::INFO <<
"RecMdcKalTrackCol registered successfully!" <<endreq;
5886 DataObject *aRecKalSegEvent;
5887 eventSvc()->findObject(
"/Event/Recon/RecMdcKalHelixSegCol", aRecKalSegEvent);
5888 if(aRecKalSegEvent!=NULL) {
5890 segsc = eventSvc()->unregisterObject(
"/Event/Recon/RecMdcKalHelixSegCol");
5891 if(segsc != StatusCode::SUCCESS) {
5892 log << MSG::FATAL <<
"Could not unregister RecMdcKalHelixSegCol collection" << endreq;
5897 segsc = eventSvc()->registerObject(
"/Event/Recon/RecMdcKalHelixSegCol", segcol);
5898 if( segsc.isFailure() ) {
5899 log << MSG::FATAL <<
"Could not register RecMdcKalHelixSeg" << endreq;
5902 log << MSG::INFO <<
"RecMdcKalHelixSegCol registered successfully!" <<endreq;
5905 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.);
5906 double r1(0.),r2(0.),kap1(999.),kap2(999.),tanl1(0.),tanl2(0.);
5909 SmartDataPtr<RecMdcKalTrackCol> kaltrkCol(eventSvc(),
"/Event/Recon/RecMdcKalTrackCol");
5911 log << MSG::FATAL <<
"Could not find RecMdcKalTrackCol" << endreq;
5914 log << MSG::INFO <<
"Begin to check RecMdcKalTrackCol"<<endreq;
5915 RecMdcKalTrackCol::iterator iter_trk = kaltrkCol->begin();
5916 for(
int jj=1; iter_trk != kaltrkCol->end(); iter_trk++,jj++) {
5917 log << MSG::DEBUG <<
"retrieved MDC Kalmantrack:"
5918 <<
"Track Id: " << (*iter_trk)->getTrackId()
5919 <<
" Mass of the fit: "<< (*iter_trk)->getMass(2)<< endreq
5920 <<
" Length of the track: "<< (*iter_trk)->getLength(2)
5921 <<
" Tof of the track: "<< (*iter_trk)->getTof(2) << endreq
5922 <<
" Chisq of the fit: "<< (*iter_trk)->getChisq(0,2)
5923 <<
" "<< (*iter_trk)->getChisq(1,2) << endreq
5924 <<
"Ndf of the fit: "<< (*iter_trk)->getNdf(0,1)
5925 <<
" "<< (*iter_trk)->getNdf(1,2) << endreq
5926 <<
"Kappa " << (*iter_trk)->getZHelix()[2]
5927 <<
"zhelixmu "<<(*iter_trk)->getZHelixMu()
5932 std::cout<<
"the size of gothelixsegs ..."<<gothelixsegs.size()<<std::endl;
5935 HelixSegRefVec::iterator it_gothelixseg = gothelixsegs.begin();
5936 for( ; it_gothelixseg != gothelixsegs.end(); it_gothelixseg++) {
5938 std::cout<<
"the layerId of this helixseg is ..."<<(*it_gothelixseg)->getLayerId()<<std::endl;
5939 std::cout<<
"the residual of this helixseg exclude the meas hit"<<(*it_gothelixseg)->getResExcl()<<std::endl;
5940 std::cout<<
"the residual of this helixseg include the meas hit"<<(*it_gothelixseg)->getResIncl()<<std::endl;
5941 std::cout<<
"the track id of the helixseg is ..."<<(*it_gothelixseg)->getTrackId() <<std::endl;
5942 std::cout<<
"the tof of the helixseg is ..."<<(*it_gothelixseg)->getTof()<<std::endl;
5943 std::cout<<
"the Zhit of the helixseg is ..."<<(*it_gothelixseg)->getZhit()<<std::endl;
5946 for(
int i = 0; i<43; i++) {
5947 log << MSG::DEBUG <<
"retrieved pathl["<<i<<
"]= "
5948 << (*iter_trk)->getPathl(i) <<endreq;
5952 m_trackid = (*iter_trk)->getTrackId();
5953 for(
int jj =0, iii=0; jj<5; jj++) {
5954 m_length[jj] = (*iter_trk)->getLength(jj);
5955 m_tof[jj] = (*iter_trk)->getTof(jj);
5956 m_nhits[jj] = (*iter_trk)->getNhits(jj);
5957 m_zhelix[jj] = (*iter_trk)->getZHelix()[jj];
5958 m_zhelixe[jj] = (*iter_trk)->getZHelixE()[jj];
5959 m_zhelixmu[jj] = (*iter_trk)->getZHelixMu()[jj];
5960 m_zhelixk[jj] = (*iter_trk)->getZHelixK()[jj];
5961 m_zhelixp[jj] = (*iter_trk)->getZHelixP()[jj];
5962 m_fhelix[jj] = (*iter_trk)->getFHelix()[jj];
5963 m_fhelixe[jj] = (*iter_trk)->getFHelixE()[jj];
5964 m_fhelixmu[jj] = (*iter_trk)->getFHelixMu()[jj];
5965 m_fhelixk[jj] = (*iter_trk)->getFHelixK()[jj];
5966 m_fhelixp[jj] = (*iter_trk)->getFHelixP()[jj];
5967 m_lhelix[jj] = (*iter_trk)->getLHelix()[jj];
5968 m_lhelixe[jj] = (*iter_trk)->getLHelixE()[jj];
5969 m_lhelixmu[jj] = (*iter_trk)->getLHelixMu()[jj];
5970 m_lhelixk[jj] = (*iter_trk)->getLHelixK()[jj];
5971 m_lhelixp[jj] = (*iter_trk)->getLHelixP()[jj];
5973 for(
int kk=0; kk<=jj; kk++,iii++) {
5974 m_zerror[iii] = (*iter_trk)->getZError()[jj][kk];
5975 m_zerrore[iii] = (*iter_trk)->getZErrorE()[jj][kk];
5976 m_zerrormu[iii] = (*iter_trk)->getZErrorMu()[jj][kk];
5977 m_zerrork[iii] = (*iter_trk)->getZErrorK()[jj][kk];
5978 m_zerrorp[iii] = (*iter_trk)->getZErrorP()[jj][kk];
5979 m_ferror[iii] = (*iter_trk)->getFError()[jj][kk];
5980 m_ferrore[iii] = (*iter_trk)->getFErrorE()[jj][kk];
5981 m_ferrormu[iii] = (*iter_trk)->getFErrorMu()[jj][kk];
5982 m_ferrork[iii] = (*iter_trk)->getFErrorK()[jj][kk];
5983 m_ferrorp[iii] = (*iter_trk)->getFErrorP()[jj][kk];
5984 m_lerror[iii] = (*iter_trk)->getLError()[jj][kk];
5985 m_lerrore[iii] = (*iter_trk)->getLErrorE()[jj][kk];
5986 m_lerrormu[iii] = (*iter_trk)->getLErrorMu()[jj][kk];
5987 m_lerrork[iii] = (*iter_trk)->getLErrorK()[jj][kk];
5988 m_lerrorp[iii] = (*iter_trk)->getLErrorP()[jj][kk];
5995 m_chisq[0][0] = (*iter_trk)->getChisq(0,0);
5996 m_chisq[1][0] = (*iter_trk)->getChisq(0,1);
5997 m_chisq[2][0] = (*iter_trk)->getChisq(0,2);
5998 m_chisq[3][0] = (*iter_trk)->getChisq(0,3);
5999 m_chisq[4][0] = (*iter_trk)->getChisq(0,4);
6000 m_chisq[0][1] = (*iter_trk)->getChisq(1,0);
6001 m_chisq[1][1] = (*iter_trk)->getChisq(1,1);
6002 m_chisq[2][1] = (*iter_trk)->getChisq(1,2);
6003 m_chisq[3][1] = (*iter_trk)->getChisq(1,3);
6004 m_chisq[4][1] = (*iter_trk)->getChisq(1,4);
6006 m_ndf[0][0] = (*iter_trk)->getNdf(0,0);
6007 m_ndf[1][0] = (*iter_trk)->getNdf(0,1);
6008 m_ndf[2][0] = (*iter_trk)->getNdf(0,2);
6009 m_ndf[3][0] = (*iter_trk)->getNdf(0,3);
6010 m_ndf[4][0] = (*iter_trk)->getNdf(0,4);
6011 m_ndf[0][1] = (*iter_trk)->getNdf(1,0);
6012 m_ndf[1][1] = (*iter_trk)->getNdf(1,1);
6013 m_ndf[2][1] = (*iter_trk)->getNdf(1,2);
6014 m_ndf[3][1] = (*iter_trk)->getNdf(1,3);
6015 m_ndf[4][1] = (*iter_trk)->getNdf(1,4);
6017 m_stat[0][0] = (*iter_trk)->getStat(0,0);
6018 m_stat[1][0] = (*iter_trk)->getStat(0,1);
6019 m_stat[2][0] = (*iter_trk)->getStat(0,2);
6020 m_stat[3][0] = (*iter_trk)->getStat(0,3);
6021 m_stat[4][0] = (*iter_trk)->getStat(0,4);
6022 m_stat[0][1] = (*iter_trk)->getStat(1,0);
6023 m_stat[1][1] = (*iter_trk)->getStat(1,1);
6024 m_stat[2][1] = (*iter_trk)->getStat(1,2);
6025 m_stat[3][1] = (*iter_trk)->getStat(1,3);
6026 m_stat[4][1] = (*iter_trk)->getStat(1,4);
6028 m_fptot = sqrt(1+pow(m_fhelix[4],2))/m_fhelix[2];
6029 m_fptote = sqrt(1+pow(m_fhelixe[4],2))/m_fhelixe[2];
6030 m_fptotmu = sqrt(1+pow(m_fhelixmu[4],2))/m_fhelixmu[2];
6031 m_fptotk = sqrt(1+pow(m_fhelixk[4],2))/m_fhelixk[2];
6032 m_fptotp = sqrt(1+pow(m_fhelixp[4],2))/m_fhelixp[2];
6034 m_zpt = 1/m_zhelix[2];
6035 m_zpte = 1/m_zhelixe[2];
6036 m_zptmu = 1/m_zhelixmu[2];
6037 m_zptk = 1/m_zhelixk[2];
6038 m_zptp = 1/m_zhelixp[2];
6040 m_fpt = 1/m_fhelix[2];
6041 m_fpte = 1/m_fhelixe[2];
6042 m_fptmu = 1/m_fhelixmu[2];
6043 m_fptk = 1/m_fhelixk[2];
6044 m_fptp = 1/m_fhelixp[2];
6046 m_lpt = 1/m_lhelix[2];
6047 m_lpte = 1/m_lhelixe[2];
6048 m_lptmu = 1/m_lhelixmu[2];
6049 m_lptk = 1/m_lhelixk[2];
6050 m_lptp = 1/m_lhelixp[2];
6052 m_lptot = sqrt(1+pow(m_lhelix[4],2))/m_lhelix[2];
6053 m_lptote = sqrt(1+pow(m_lhelixe[4],2))/m_lhelixe[2];
6054 m_lptotmu = sqrt(1+pow(m_lhelixmu[4],2))/m_lhelixmu[2];
6055 m_lptotk = sqrt(1+pow(m_lhelixk[4],2))/m_lhelixk[2];
6056 m_lptotp = sqrt(1+pow(m_lhelixp[4],2))/m_lhelixp[2];
6058 m_zptot = sqrt(1+pow(m_zhelix[4],2))/m_zhelix[2];
6059 m_zptote = sqrt(1+pow(m_zhelixe[4],2))/m_zhelixe[2];
6060 m_zptotmu = sqrt(1+pow(m_zhelixmu[4],2))/m_zhelixmu[2];
6061 m_zptotk = sqrt(1+pow(m_zhelixk[4],2))/m_zhelixk[2];
6062 m_zptotp = sqrt(1+pow(m_zhelixp[4],2))/m_zhelixp[2];
6064 m_zsigp = sqrt(pow((m_zptot/m_zhelix[2]),2)*m_zerror[5]+
6065 pow((m_zhelix[4]/m_zptot),2)*pow((1/m_zhelix[2]),4)*m_zerror[14]-
6066 2*m_zhelix[4]*m_zerror[12]*pow((1/m_zhelix[2]),3));
6067 m_zsigpe = sqrt(pow((m_zptote/m_zhelixe[2]),2)*m_zerrore[5]+
6068 pow((m_zhelixe[4]/m_zptote),2)*pow((1/m_zhelixe[2]),4)*m_zerrore[14]-
6069 2*m_zhelixe[4]*m_zerrore[12]*pow((1/m_zhelixe[2]),3));
6070 m_zsigpmu = sqrt(pow((m_zptotmu/m_zhelixmu[2]),2)*m_zerrormu[5]+
6071 pow((m_zhelixmu[4]/m_zptotmu),2)*pow((1/m_zhelixmu[2]),4)*m_zerrormu[14]-
6072 2*m_zhelixmu[4]*m_zerrormu[12]*pow((1/m_zhelixmu[2]),3));
6073 m_zsigpk = sqrt(pow((m_zptotk/m_zhelixk[2]),2)*m_zerrork[5]+
6074 pow((m_zhelixk[4]/m_zptotk),2)*pow((1/m_zhelixk[2]),4)*m_zerrork[14]-
6075 2*m_zhelixk[4]*m_zerrork[12]*pow((1/m_zhelixk[2]),3));
6076 m_zsigpp = sqrt(pow((m_zptotp/m_zhelixp[2]),2)*m_zerrorp[5]+
6077 pow((m_zhelixp[4]/m_zptotp),2)*pow((1/m_zhelixp[2]),4)*m_zerrorp[14]-
6078 2*m_zhelixp[4]*m_zerrorp[12]*pow((1/m_zhelixp[2]),3));
6081 StatusCode sc1 = m_nt1->write();
6082 if( sc1.isFailure() ) cout<<
"Ntuple1 filling failed!"<<endl;
6087 phi1 = (*iter_trk)->getFFi0();
6088 r1 = (*iter_trk)->getFDr();
6089 z1 = (*iter_trk)->getFDz();
6090 kap1 = (*iter_trk)->getFCpa();
6091 tanl1 = (*iter_trk)->getFTanl();
6094 p1 = sqrt(1+tanl1*tanl1)/kap1;
6095 the1 =
M_PI/2-atan(tanl1);
6096 }
else if(jj == 2) {
6097 phi2 = (*iter_trk)->getFFi0();
6098 r2 = (*iter_trk)->getFDr();
6099 z2 = (*iter_trk)->getFDz();
6100 kap2 = (*iter_trk)->getFCpa();
6101 tanl2 = (*iter_trk)->getFTanl();
6104 p2 = sqrt(1+tanl2*tanl2)/kap1;
6105 the2 =
M_PI/2-atan(tanl2);
6113 m_delthe = the1 + the2;
6116 StatusCode sc2 = m_nt2->write();
6117 if( sc2.isFailure() ) cout<<
"Ntuple2 filling failed!"<<endl;
6125 cout <<
"Kalfitting finished " << std::endl;
7387 std::map<double,double> esti1_r;
7388 std::map<double,double> esti2_r;
7389 std::map<double,double> esti1_phi;
7390 std::map<double,double> esti2_phi;
7391 std::map<double,double> esti1_z;
7392 std::map<double,double> esti2_z;
7393 std::map<double,double> diff1_phi;
7394 std::map<double,double> diff1_z;
7395 std::map<double,double> diff2_phi;
7396 std::map<double,double> diff2_z;
7397 std::map<double,HepPoint3D> posnew;
7398 std::map<double,HepSymMatrix> eanew;
7399 std::map<double,HepVector> anew;
7400 std::map<double,double> meas_r;
7401 std::map<double,double> meas_phi;
7402 std::map<double,double> meas_z;
7403 std::map<double,double> meas_phierr;
7404 std::map<double,double> esti_phierr;
7405 std::map<double,double> meas_zerr;
7406 std::map<double,double> esti_zerr;
7408 if(
debug_==4){ cout <<
"wall size" << _BesKalmanFitWalls.size() << endl;}
7409 for(
int iLayer=2; iLayer>=0; iLayer--)
7427 meas_phierr.clear();
7429 esti_phierr.clear();
7432 int nHits=myGemHitCol[iLayer].size();
7435 if (nHits <= 0)
continue;
7438 vector<KalFitGemHit>::iterator
iter;
7442 innerwall(track,hypo,way,64+iLayer*(-32),88+iLayer*(-32));
7445 for(
iter=myGemHitCol[iLayer].begin();
iter!=myGemHitCol[iLayer].end();
iter++ )
7448 track.
set(track_start.
pivot(),track_start.
a(),track_start.
Ea());
7450 double rGem=(*iter).getR();
7451 double phiGem=(*iter).getPhi();
7452 double zGem=(*iter).getZ();
7453 HepVector v_measu=(*iter).getVecPos();
7466 double x0=posEstiAtGem.x();
7467 double y0=posEstiAtGem.y();
7475 if(
muls_) track.
ms(pathInGem, _BesKalmanFitMaterials[26+32*iLayer], way);
7476 if(
loss_) track.
eloss(pathInGem, _BesKalmanFitMaterials[26+32*iLayer], way);
7479 double rGemEsti=posEstiAtGem.perp();
7480 double phiGemEsti=posEstiAtGem.phi();
7481 double zGemEsti=posEstiAtGem.z();
7482 HepVector v_estim(2,0);
7483 v_estim(1)=phiGemEsti;
7484 v_estim(2)=zGemEsti;
7491 const HepSymMatrix& Ea = track.
Ea();
7492 HepVector v_a = track.
a();
7495 double kappa=v_a(3);
7500 HepMatrix
H(2, 5, 0);
7503 H(1,1)=-y0*
cos(phi0)/(y0*y0+x0*x0)+x0*
sin(phi0)/(x0*x0+y0*y0);
7506 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)));
7507 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));
7508 H(2,3)=
m_alpha/(kappa*kappa)*tl*dPhi;
7512 HepMatrix H_T=
H.T();
7517 HepSymMatrix V=(*iter).getErrMat();
7520 HepMatrix HEaH_T=HEa*H_T;
7522 HepMatrix Vinv=(V+HEaH_T).inverse(ierr);
7523 if(ierr!=0) cout<<
"error in HEaH_T.inverse()!"<<endl;
7525 HepMatrix K=Ea*H_T*Vinv;
7527 HepSymMatrix EaNew(5,0);
7528 EaNew.assign(Ea-K*
H*Ea);
7530 HepVector v_diff=v_measu-v_estim;
7532 double delPhi=v_diff(1);
7533 if(fabs(v_diff(1))>6.28) v_measu(1)=-1*v_measu(1);
7540 HepVector aNew=v_a+K*v_diff;
7546 HepPoint3D posEstiAtGem_new = track_test.
x(dPhiToGem_New);
7547 double phiGemEsti_new=posEstiAtGem_new.phi();
7548 double zGemEsti_new=posEstiAtGem_new.z();
7549 double x0_new = posEstiAtGem_new.x();
7550 double y0_new = posEstiAtGem_new.y();
7551 HepVector v_estim_new(2,0);
7552 v_estim_new(1)=phiGemEsti_new;
7553 v_estim_new(2)=zGemEsti_new;
7557 v_diff=v_measu-v_estim_new;
7565 track_test.
pivot_numf(posEstiAtGem_new, pathInGem);
7566 HepVector v_a_new = track_test.
a();
7568 HepSymMatrix Ea_new = track_test.
Ea();
7570 double phi0_new = v_a_new(2);
7571 double kappa_new=v_a_new(3);
7572 HepMatrix H_new(2, 5, 0);
7581 HepMatrix H_new_T=H_new.T();
7588 HepSymMatrix R(2,0);
7590 R.assign(V-
H*EaNew*H_T);
7598 HepVector v_dChi2=v_diff.T()*R.inverse(ierr)*v_diff;
7599 if(ierr!=0) cout<<
"error in R.inverse()!"<<endl;
7602 double dChi2=v_dChi2(1);
7606 cout<<
"pivot: "<<posEstiAtGem<<endl;
7607 cout<<
"new pivot: "<<posEstiAtGem<<endl;
7608 cout<<
"a: "<<v_a<<endl;
7609 cout<<
"new a: "<<aNew<<endl;
7610 cout<<
"Ea: "<<Ea<<endl;
7611 cout<<
"new Ea: "<<EaNew<<endl;
7612 cout<<
"dchi2= "<<dChi2<<endl;
7647 esti1_r[dChi2]=rGemEsti;
7648 esti1_phi[dChi2]=phiGemEsti;
7649 esti1_z[dChi2]=zGemEsti;
7650 esti2_r[dChi2]=rGem;
7651 esti2_phi[dChi2]=phiGemEsti_new;
7652 esti2_z[dChi2]=zGemEsti_new;
7653 diff1_phi[dChi2]=phiGem-phiGemEsti;
7654 diff1_z[dChi2]=zGem-zGemEsti;
7655 diff2_phi[dChi2]=phiGem-phiGemEsti_new;
7656 diff2_z[dChi2]=zGem-zGemEsti_new;
7657 posnew[dChi2]=posEstiAtGem;
7661 meas_phi[dChi2]=phiGem;
7664 meas_phierr[dChi2]=sqrt(V[0][0]);
7666 meas_zerr[dChi2]=sqrt(V[1][1]);
7668 esti_phierr[dChi2]=sqrt((
H*(track_start.
Ea())*H_T)[0][0]);
7670 esti_zerr[dChi2]=sqrt((
H*(track_start.
Ea())*H_T)[1][1]);
7679 track.
set((posnew.begin())->second,(anew.begin())->second,(eanew.begin()->second));
7683 HepPoint3D new_posEstiAtGem = track.
x(new_dPhiToGem);
7684 double new_pathInGem;
7687 track.
pivot_numf(new_posEstiAtGem, new_pathInGem);
7692 if(
muls_) track.
ms(new_pathInGem, _BesKalmanFitMaterials[26+32*iLayer], way);
7693 if(
loss_) track.
eloss(new_pathInGem, _BesKalmanFitMaterials[26+32*iLayer], way);
7697 innerwall(track,hypo,way,90+(-32)*iLayer,95+(-32)*iLayer);}
7699 innerwall(track,hypo,way,90+(-32)*iLayer,94+(-32)*iLayer);}
7707 m_esti1_r[ii]=(esti1_r.begin())->second;
7709 m_esti1_phi[ii]=(esti1_phi.begin())->second;
7710 m_esti1_z[ii]=(esti1_z.begin())->second;
7711 m_esti2_r[ii]=(esti2_r.begin())->second;
7712 m_esti2_phi[ii]=(esti2_phi.begin())->second;
7713 m_esti2_z[ii]=(esti2_z.begin())->second;
7714 m_diff1_phi[ii]=(diff1_phi.begin())->second;
7715 m_diff1_z[ii]=(diff1_z.begin())->second;
7716 m_diff2_phi[ii]=(diff2_phi.begin())->second;
7717 m_diff2_z[ii]=(diff2_z.begin())->second;
7718 m_Gchi2[ii]=(esti1_r.begin())->first;
7719 m_meas_r[ii]=(meas_r.begin())->second;
7720 m_meas_phi[ii]=(meas_phi.begin())->second;
7722 m_meas_z[ii]=(meas_z.begin())->second;
7723 m_meas_phierr[ii]=(meas_phierr.begin())->second;
7724 m_meas_zerr[ii]=(meas_zerr.begin())->second;
7725 m_esti_phierr[ii]=(esti_phierr.begin())->second;
7726 m_esti_zerr[ii]=(esti_zerr.begin())->second;
7727 m_GemLayer[ii]=iLayer;
7739 if(m_meas_phi[0]>0 && m_meas_phi[0]<1.5707963) m_qua=0;
7740 if(m_meas_phi[0]>1.5707963 && m_meas_phi[0]<3.1415926) m_qua=1;
7741 if(m_meas_phi[0]>-3.1415926&& m_meas_phi[0]<-1.5707963) m_qua=2;
7742 if(m_meas_phi[0]>-1.5707963&& m_meas_phi[0]<0) m_qua=3;
7765 meas_phierr.clear();
7767 esti_phierr.clear();
7773 SmartDataPtr<Event::EventHeader> evtHeader(eventSvc(),
"/Event/EventHeader");
7775 m_evt3 = evtHeader->eventNumber();
7787 SmartDataPtr<Event::EventHeader> eventHeader(eventSvc(),
"/Event/EventHeader");
7789 cout <<
"Could not find Event Header" << endl;
7791 int eventNo = eventHeader->eventNumber();
7792 int runNo = eventHeader->runNumber();
7793 SmartDataPtr<RecMdcTrackCol> recMdcTrack(eventSvc(),
"/Event/Recon/RecMdcTrackCol");
7795 cout <<
"Could not find RecMdcTrackCol" << endl;
7797 RecMdcTrackCol::iterator itrk = recMdcTrack->begin();
7799 for(; itrk != recMdcTrack->end(); itrk++){
7800 if((*itrk)->trackId() == track.
trasan_id() && (*itrk)->getNcluster() > 0) {
7801 clusterRefVec = (*itrk)->getVecClusters();
7809 if(
debug_==4){ cout<<
"wall size"<<_BesKalmanFitWalls.size()<<endl;}
7812 ClusterRefVec::iterator itCluster;
7813 ClusterRefVec::iterator itStartFlag;
7814 ClusterRefVec::iterator itStopFlag;
7818 double R_o_cgem=_BesKalmanFitWalls[0].radius();
7820 itStartFlag = clusterRefVec.begin();
7821 itStopFlag = clusterRefVec.end();
7826 itStartFlag = clusterRefVec.end()-1;
7827 itStopFlag = clusterRefVec.begin()-1;
7837 for(itCluster = itStartFlag; itCluster !=itStopFlag; itCluster = itCluster + step){
7838 int layer = (*itCluster)->getlayerid();
7841 innerwall(track,hypo,way,lastR,dmR,index);
7851 m_dchi2_hit_e[m_hit_no] = dchi2;
7852 m_residest_hit_e[m_hit_no] = 0;
7853 m_residnew_hit_e[m_hit_no] = 0;
7854 m_layer_hit_e[m_hit_no] = layer;
7856 m_anal_dr_hit_e[m_hit_no] = worktemp.
a()[0];
7857 m_anal_phi0_hit_e[m_hit_no] = worktemp.
a()[1];
7858 m_anal_kappa_hit_e[m_hit_no] = worktemp.
a()[2];
7859 m_anal_dz_hit_e[m_hit_no] = worktemp.
a()[3];
7860 m_anal_tanl_hit_e[m_hit_no] = worktemp.
a()[4];
7861 m_anal_ea_dr_hit_e[m_hit_no] = worktemp.
Ea()[0][0];
7862 m_anal_ea_phi0_hit_e[m_hit_no] = worktemp.
Ea()[1][1];
7863 m_anal_ea_kappa_hit_e[m_hit_no] = worktemp.
Ea()[2][2];
7864 m_anal_ea_dz_hit_e[m_hit_no] = worktemp.
Ea()[3][3];
7865 m_anal_ea_tanl_hit_e[m_hit_no] = worktemp.
Ea()[4][4];
7870 m_dchi2_hit_mu[m_hit_no] = dchi2;
7871 m_residest_hit_mu[m_hit_no] = 0;
7872 m_residnew_hit_mu[m_hit_no] = 0;
7873 m_layer_hit_mu[m_hit_no] = layer;
7875 m_anal_dr_hit_mu[m_hit_no] = worktemp.
a()[0];
7876 m_anal_phi0_hit_mu[m_hit_no] = worktemp.
a()[1];
7877 m_anal_kappa_hit_mu[m_hit_no] = worktemp.
a()[2];
7878 m_anal_dz_hit_mu[m_hit_no] = worktemp.
a()[3];
7879 m_anal_tanl_hit_mu[m_hit_no] = worktemp.
a()[4];
7880 m_anal_ea_dr_hit_mu[m_hit_no] = worktemp.
Ea()[0][0];
7881 m_anal_ea_phi0_hit_mu[m_hit_no] = worktemp.
Ea()[1][1];
7882 m_anal_ea_kappa_hit_mu[m_hit_no] = worktemp.
Ea()[2][2];
7883 m_anal_ea_dz_hit_mu[m_hit_no] = worktemp.
Ea()[3][3];
7884 m_anal_ea_tanl_hit_mu[m_hit_no] = worktemp.
Ea()[4][4];
7890 m_dchi2_hit_pi[m_hit_no] = dchi2;
7891 m_residest_hit_pi[m_hit_no] = 0;
7892 m_residnew_hit_pi[m_hit_no] = 0;
7893 m_layer_hit_pi[m_hit_no] = layer;
7895 m_anal_dr_hit_pi[m_hit_no] = worktemp.
a()[0];
7896 m_anal_phi0_hit_pi[m_hit_no] = worktemp.
a()[1];
7897 m_anal_kappa_hit_pi[m_hit_no] = worktemp.
a()[2];
7898 m_anal_dz_hit_pi[m_hit_no] = worktemp.
a()[3];
7899 m_anal_tanl_hit_pi[m_hit_no] = worktemp.
a()[4];
7900 m_anal_ea_dr_hit_pi[m_hit_no] = worktemp.
Ea()[0][0];
7901 m_anal_ea_phi0_hit_pi[m_hit_no] = worktemp.
Ea()[1][1];
7902 m_anal_ea_kappa_hit_pi[m_hit_no] = worktemp.
Ea()[2][2];
7903 m_anal_ea_dz_hit_pi[m_hit_no] = worktemp.
Ea()[3][3];
7904 m_anal_ea_tanl_hit_pi[m_hit_no] = worktemp.
Ea()[4][4];
7909 m_dchi2_hit_k[m_hit_no] = dchi2;
7910 m_residest_hit_k[m_hit_no] = 0;
7911 m_residnew_hit_k[m_hit_no] = 0;
7912 m_layer_hit_k[m_hit_no] = layer;
7914 m_anal_dr_hit_k[m_hit_no] = worktemp.
a()[0];
7915 m_anal_phi0_hit_k[m_hit_no] = worktemp.
a()[1];
7916 m_anal_kappa_hit_k[m_hit_no] = worktemp.
a()[2];
7917 m_anal_dz_hit_k[m_hit_no] = worktemp.
a()[3];
7918 m_anal_tanl_hit_k[m_hit_no] = worktemp.
a()[4];
7919 m_anal_ea_dr_hit_k[m_hit_no] = worktemp.
Ea()[0][0];
7920 m_anal_ea_phi0_hit_k[m_hit_no] = worktemp.
Ea()[1][1];
7921 m_anal_ea_kappa_hit_k[m_hit_no] = worktemp.
Ea()[2][2];
7922 m_anal_ea_dz_hit_k[m_hit_no] = worktemp.
Ea()[3][3];
7923 m_anal_ea_tanl_hit_k[m_hit_no] = worktemp.
Ea()[4][4];
7928 m_dchi2_hit_p[m_hit_no] = dchi2;
7929 m_residest_hit_p[m_hit_no] = 0;
7930 m_residnew_hit_p[m_hit_no] = 0;
7931 m_layer_hit_p[m_hit_no] = layer;
7933 m_anal_dr_hit_p[m_hit_no] = worktemp.
a()[0];
7934 m_anal_phi0_hit_p[m_hit_no] = worktemp.
a()[1];
7935 m_anal_kappa_hit_p[m_hit_no] = worktemp.
a()[2];
7936 m_anal_dz_hit_p[m_hit_no] = worktemp.
a()[3];
7937 m_anal_tanl_hit_p[m_hit_no] = worktemp.
a()[4];
7938 m_anal_ea_dr_hit_p[m_hit_no] = worktemp.
Ea()[0][0];
7939 m_anal_ea_phi0_hit_p[m_hit_no] = worktemp.
Ea()[1][1];
7940 m_anal_ea_kappa_hit_p[m_hit_no] = worktemp.
Ea()[2][2];
7941 m_anal_ea_dz_hit_p[m_hit_no] = worktemp.
Ea()[3][3];
7942 m_anal_ea_tanl_hit_p[m_hit_no] = worktemp.
Ea()[4][4];
7950 innerwall(track,hypo,way,lastR,R_o_cgem,index);
7954 innerwall(track,hypo,way,lastR,0,index);