494 auto sink = FairRunAna::Instance()->GetSink();
495 assert(sink->GetSinkType() == kFILESINK);
496 auto rootFileSink =
static_cast<FairRootFileSink*
>(sink);
497 fHist = rootFileSink->GetRootFile();
499 TString hname =
fHist->GetName();
500 hname.Insert(hname.Length() - 5,
".HadAna");
501 fHist =
new TFile(hname,
"recreate");
502 LOG(info) <<
"CreateHistograms in " <<
fHist->GetName();
504 fHist->ReOpen(
"Update");
513 Float_t ptmmax = 2.5;
520 new TH2F(
"ptm_rap_gen_pip",
"MCTrack-gen pi-plus; y; p_{T}/m", ptm_nbx, ymin, ymax, ptm_nby, 0., ptmmax);
522 new TH2F(
"ptm_rap_gen_pim",
"MCTrack-gen pi-minus;y; p_{T}/m", ptm_nbx, ymin, ymax, ptm_nby, 0., ptmmax);
524 new TH2F(
"ptm_rap_gen_kp",
"MCTrack-gen k-plus; y; p_{T}/m", ptm_nbx, ymin, ymax, ptm_nby, 0., ptmmax);
526 new TH2F(
"ptm_rap_gen_km",
"MCTrack-gen k-minus; y; p_{T}/m", ptm_nbx, ymin, ymax, ptm_nby, 0., ptmmax);
528 new TH2F(
"ptm_rap_gen_p",
"MCTrack-gen proton; y; p_{T}/m", ptm_nbx, ymin, ymax, ptm_nby, 0., ptmmax);
530 new TH2F(
"ptm_rap_gen_pbar",
"MCTrack-gen antiproton;y;p_{T}/m", ptm_nbx, ymin, ymax, ptm_nby, 0., ptmmax);
532 new TH2F(
"ptm_rap_gen_d",
"MCTrack-gen deuteron;y; p_{T}/m", ptm_nbx, ymin, ymax, ptm_nby, 0., ptmmax);
534 new TH2F(
"ptm_rap_gen_t",
"MCTrack-gen triton; y; p_{T}/m", ptm_nbx, ymin, ymax, ptm_nby, 0., ptmmax);
536 new TH2F(
"ptm_rap_gen_h",
"MCTrack-gen 3he; y; p_{T}/m", ptm_nbx, ymin, ymax, ptm_nby, 0., ptmmax);
538 new TH2F(
"ptm_rap_gen_a",
"MCTrack-gen alpha; y; p_{T}/m", ptm_nbx, ymin, ymax, ptm_nby, 0., ptmmax);
540 new TH2F(
"ptm_rap_gen_imf",
"MCTrack-gen imf; y; p_{T}/m", ptm_nbx, ymin, ymax, ptm_nby, 0., ptmmax);
542 Float_t v1_nbx = 20.;
543 Float_t v1_nby = 20.;
547 new TH2F(
"v1_rap_gen_pip",
"MCTrack-gen pi-plus; y; v_{1}", v1_nbx, -yvmax, yvmax, v1_nby, -1., 1.);
549 new TH2F(
"v1_rap_gen_pim",
"MCTrack-gen pi-minus;y; v_{1}", v1_nbx, -yvmax, yvmax, v1_nby, -1., 1.);
550 fa_v1_rap_gen_kp =
new TH2F(
"v1_rap_gen_kp",
"MCTrack-gen k-plus; y; v_{1}", v1_nbx, -yvmax, yvmax, v1_nby, -1., 1.);
551 fa_v1_rap_gen_km =
new TH2F(
"v1_rap_gen_km",
"MCTrack-gen k-minus; y; v_{1}", v1_nbx, -yvmax, yvmax, v1_nby, -1., 1.);
552 fa_v1_rap_gen_p =
new TH2F(
"v1_rap_gen_p",
"MCTrack-gen proton; y; v_{1}", v1_nbx, -yvmax, yvmax, v1_nby, -1., 1.);
554 new TH2F(
"v1_rap_gen_pbar",
"MCTrack-gen antiproton;y; v_{1}", v1_nbx, -yvmax, yvmax, v1_nby, -1., 1.);
555 fa_v1_rap_gen_d =
new TH2F(
"v1_rap_gen_d",
"MCTrack-gen deuteron;y; v_{1}", v1_nbx, -yvmax, yvmax, v1_nby, -1., 1.);
556 fa_v1_rap_gen_t =
new TH2F(
"v1_rap_gen_t",
"MCTrack-gen triton; y; v_{1}", v1_nbx, -yvmax, yvmax, v1_nby, -1., 1.);
557 fa_v1_rap_gen_h =
new TH2F(
"v1_rap_gen_h",
"MCTrack-gen 3he; y; v_{1}", v1_nbx, -yvmax, yvmax, v1_nby, -1., 1.);
558 fa_v1_rap_gen_a =
new TH2F(
"v1_rap_gen_a",
"MCTrack-gen alpha; y; v_{1}", v1_nbx, -yvmax, yvmax, v1_nby, -1., 1.);
560 new TH2F(
"v1_rap_gen_imf",
"MCTrack-gen imf; y; v_{1}", v1_nbx, -yvmax, yvmax, v1_nby, -1., 1.);
563 new TH2F(
"v2_rap_gen_pip",
"MCTrack-gen pi-plus; y; v_{2}", v1_nbx, -yvmax, yvmax, v1_nby, -1., 1.);
565 new TH2F(
"v2_rap_gen_pim",
"MCTrack-gen pi-minus;y; v_{2}", v1_nbx, -yvmax, yvmax, v1_nby, -1., 1.);
566 fa_v2_rap_gen_kp =
new TH2F(
"v2_rap_gen_kp",
"MCTrack-gen k-plus; y; v_{2}", v1_nbx, -yvmax, yvmax, v1_nby, -1., 1.);
567 fa_v2_rap_gen_km =
new TH2F(
"v2_rap_gen_km",
"MCTrack-gen k-minus; y; v_{2}", v1_nbx, -yvmax, yvmax, v1_nby, -1., 1.);
568 fa_v2_rap_gen_p =
new TH2F(
"v2_rap_gen_p",
"MCTrack-gen proton; y; v_{2}", v1_nbx, -yvmax, yvmax, v1_nby, -1., 1.);
570 new TH2F(
"v2_rap_gen_pbar",
"MCTrack-gen antiproton;y; v_{2}", v1_nbx, -yvmax, yvmax, v1_nby, -1., 1.);
571 fa_v2_rap_gen_d =
new TH2F(
"v2_rap_gen_d",
"MCTrack-gen deuteron;y; v_{2}", v1_nbx, -yvmax, yvmax, v1_nby, -1., 1.);
572 fa_v2_rap_gen_t =
new TH2F(
"v2_rap_gen_t",
"MCTrack-gen triton; y; v_{2}", v1_nbx, -yvmax, yvmax, v1_nby, -1., 1.);
573 fa_v2_rap_gen_h =
new TH2F(
"v2_rap_gen_h",
"MCTrack-gen 3he; y; v_{2}", v1_nbx, -yvmax, yvmax, v1_nby, -1., 1.);
574 fa_v2_rap_gen_a =
new TH2F(
"v2_rap_gen_a",
"MCTrack-gen alpha; y; v_{2}", v1_nbx, -yvmax, yvmax, v1_nby, -1., 1.);
576 new TH2F(
"v2_rap_gen_imf",
"MCTrack-gen imf; y; v_{2}", v1_nbx, -yvmax, yvmax, v1_nby, -1., 1.);
581 new TH2F(
"ptm_rap_poi_pip",
"MCTrack-poi pi-plus; y; p_{T}/m", ptm_nbx, ymin, ymax, ptm_nby, 0., ptmmax);
583 new TH2F(
"ptm_rap_poi_pim",
"MCTrack-poi pi-minus;y; p_{T}/m", ptm_nbx, ymin, ymax, ptm_nby, 0., ptmmax);
585 new TH2F(
"ptm_rap_poi_kp",
"MCTrack-poi k-plus; y; p_{T}/m", ptm_nbx, ymin, ymax, ptm_nby, 0., ptmmax);
587 new TH2F(
"ptm_rap_poi_km",
"MCTrack-poi k-minus; y; p_{T}/m", ptm_nbx, ymin, ymax, ptm_nby, 0., ptmmax);
589 new TH2F(
"ptm_rap_poi_p",
"MCTrack-poi proton; y; p_{T}/m", ptm_nbx, ymin, ymax, ptm_nby, 0., ptmmax);
591 new TH2F(
"ptm_rap_poi_pbar",
"MCTrack-poi antiproton;y;p_{T}/m", ptm_nbx, ymin, ymax, ptm_nby, 0., ptmmax);
593 new TH2F(
"ptm_rap_poi_d",
"MCTrack-poi deuteron;y; p_{T}/m", ptm_nbx, ymin, ymax, ptm_nby, 0., ptmmax);
595 new TH2F(
"ptm_rap_poi_t",
"MCTrack-poi triton; y; p_{T}/m", ptm_nbx, ymin, ymax, ptm_nby, 0., ptmmax);
597 new TH2F(
"ptm_rap_poi_h",
"MCTrack-poi 3he; y; p_{T}/m", ptm_nbx, ymin, ymax, ptm_nby, 0., ptmmax);
599 new TH2F(
"ptm_rap_poi_a",
"MCTrack-poi alpha; y; p_{T}/m", ptm_nbx, ymin, ymax, ptm_nby, 0., ptmmax);
601 new TH2F(
"ptm_rap_poi_imf",
"MCTrack-poi imf; y; p_{T}/m", ptm_nbx, ymin, ymax, ptm_nby, 0., ptmmax);
604 new TH2F(
"v1_rap_poi_pip",
"MCTrack-poi pi-plus; y; v_{1}", v1_nbx, -yvmax, yvmax, v1_nby, -1., 1.);
606 new TH2F(
"v1_rap_poi_pim",
"MCTrack-poi pi-minus;y; v_{1}", v1_nbx, -yvmax, yvmax, v1_nby, -1., 1.);
607 fa_v1_rap_poi_kp =
new TH2F(
"v1_rap_poi_kp",
"MCTrack-poi k-plus; y; v_{1}", v1_nbx, -yvmax, yvmax, v1_nby, -1., 1.);
608 fa_v1_rap_poi_km =
new TH2F(
"v1_rap_poi_km",
"MCTrack-poi k-minus; y; v_{1}", v1_nbx, -yvmax, yvmax, v1_nby, -1., 1.);
609 fa_v1_rap_poi_p =
new TH2F(
"v1_rap_poi_p",
"MCTrack-poi proton; y; v_{1}", v1_nbx, -yvmax, yvmax, v1_nby, -1., 1.);
611 new TH2F(
"v1_rap_poi_pbar",
"MCTrack-poi antiproton;y; v_{1}", v1_nbx, -yvmax, yvmax, v1_nby, -1., 1.);
612 fa_v1_rap_poi_d =
new TH2F(
"v1_rap_poi_d",
"MCTrack-poi deuteron;y; v_{1}", v1_nbx, -yvmax, yvmax, v1_nby, -1., 1.);
613 fa_v1_rap_poi_t =
new TH2F(
"v1_rap_poi_t",
"MCTrack-poi triton; y; v_{1}", v1_nbx, -yvmax, yvmax, v1_nby, -1., 1.);
614 fa_v1_rap_poi_h =
new TH2F(
"v1_rap_poi_h",
"MCTrack-poi 3he; y; v_{1}", v1_nbx, -yvmax, yvmax, v1_nby, -1., 1.);
615 fa_v1_rap_poi_a =
new TH2F(
"v1_rap_poi_a",
"MCTrack-poi alpha; y; v_{1}", v1_nbx, -yvmax, yvmax, v1_nby, -1., 1.);
617 new TH2F(
"v1_rap_poi_imf",
"MCTrack-poi imf; y; v_{1}", v1_nbx, -yvmax, yvmax, v1_nby, -1., 1.);
620 new TH2F(
"v2_rap_poi_pip",
"MCTrack-poi pi-plus; y; v_{2}", v1_nbx, -yvmax, yvmax, v1_nby, -1., 1.);
622 new TH2F(
"v2_rap_poi_pim",
"MCTrack-poi pi-minus;y; v_{2}", v1_nbx, -yvmax, yvmax, v1_nby, -1., 1.);
623 fa_v2_rap_poi_kp =
new TH2F(
"v2_rap_poi_kp",
"MCTrack-poi k-plus; y; v_{2}", v1_nbx, -yvmax, yvmax, v1_nby, -1., 1.);
624 fa_v2_rap_poi_km =
new TH2F(
"v2_rap_poi_km",
"MCTrack-poi k-minus; y; v_{2}", v1_nbx, -yvmax, yvmax, v1_nby, -1., 1.);
625 fa_v2_rap_poi_p =
new TH2F(
"v2_rap_poi_p",
"MCTrack-poi proton; y; v_{2}", v1_nbx, -yvmax, yvmax, v1_nby, -1., 1.);
627 new TH2F(
"v2_rap_poi_pbar",
"MCTrack-poi antiproton;y; v_{2}", v1_nbx, -yvmax, yvmax, v1_nby, -1., 1.);
628 fa_v2_rap_poi_d =
new TH2F(
"v2_rap_poi_d",
"MCTrack-poi deuteron;y; v_{2}", v1_nbx, -yvmax, yvmax, v1_nby, -1., 1.);
629 fa_v2_rap_poi_t =
new TH2F(
"v2_rap_poi_t",
"MCTrack-poi triton; y; v_{2}", v1_nbx, -yvmax, yvmax, v1_nby, -1., 1.);
630 fa_v2_rap_poi_h =
new TH2F(
"v2_rap_poi_h",
"MCTrack-poi 3he; y; v_{2}", v1_nbx, -yvmax, yvmax, v1_nby, -1., 1.);
631 fa_v2_rap_poi_a =
new TH2F(
"v2_rap_poi_a",
"MCTrack-poi alpha; y; v_{2}", v1_nbx, -yvmax, yvmax, v1_nby, -1., 1.);
633 new TH2F(
"v2_rap_poi_imf",
"MCTrack-poi imf; y; v_{2}", v1_nbx, -yvmax, yvmax, v1_nby, -1., 1.);
638 new TH2F(
"ptm_rap_hit_pip",
"MCTrack-hit pi-plus; y; p_{T}/m", ptm_nbx, ymin, ymax, ptm_nby, 0., ptmmax);
640 new TH2F(
"ptm_rap_hit_pim",
"MCTrack-hit pi-minus;y; p_{T}/m", ptm_nbx, ymin, ymax, ptm_nby, 0., ptmmax);
642 new TH2F(
"ptm_rap_hit_kp",
"MCTrack-hit k-plus; y; p_{T}/m", ptm_nbx, ymin, ymax, ptm_nby, 0., ptmmax);
644 new TH2F(
"ptm_rap_hit_km",
"MCTrack-hit k-minus; y; p_{T}/m", ptm_nbx, ymin, ymax, ptm_nby, 0., ptmmax);
646 new TH2F(
"ptm_rap_hit_p",
"MCTrack-hit proton; y; p_{T}/m", ptm_nbx, ymin, ymax, ptm_nby, 0., ptmmax);
648 new TH2F(
"ptm_rap_hit_pbar",
"MCTrack-hit antiproton;y;p_{T}/m", ptm_nbx, ymin, ymax, ptm_nby, 0., ptmmax);
650 new TH2F(
"ptm_rap_hit_d",
"MCTrack-hit deuteron;y; p_{T}/m", ptm_nbx, ymin, ymax, ptm_nby, 0., ptmmax);
652 new TH2F(
"ptm_rap_hit_t",
"MCTrack-hit triton; y; p_{T}/m", ptm_nbx, ymin, ymax, ptm_nby, 0., ptmmax);
654 new TH2F(
"ptm_rap_hit_h",
"MCTrack-hit 3he; y; p_{T}/m", ptm_nbx, ymin, ymax, ptm_nby, 0., ptmmax);
656 new TH2F(
"ptm_rap_hit_a",
"MCTrack-hit alpha; y; p_{T}/m", ptm_nbx, ymin, ymax, ptm_nby, 0., ptmmax);
658 new TH2F(
"ptm_rap_hit_imf",
"MCTrack-hit imf; y; p_{T}/m", ptm_nbx, ymin, ymax, ptm_nby, 0., ptmmax);
661 new TH2F(
"v1_rap_hit_pip",
"MCTrack-hit pi-plus; y; v_{1}", v1_nbx, -yvmax, yvmax, v1_nby, -1., 1.);
663 new TH2F(
"v1_rap_hit_pim",
"MCTrack-hit pi-minus;y; v_{1}", v1_nbx, -yvmax, yvmax, v1_nby, -1., 1.);
664 fa_v1_rap_hit_kp =
new TH2F(
"v1_rap_hit_kp",
"MCTrack-hit k-plus; y; v_{1}", v1_nbx, -yvmax, yvmax, v1_nby, -1., 1.);
665 fa_v1_rap_hit_km =
new TH2F(
"v1_rap_hit_km",
"MCTrack-hit k-minus; y; v_{1}", v1_nbx, -yvmax, yvmax, v1_nby, -1., 1.);
666 fa_v1_rap_hit_p =
new TH2F(
"v1_rap_hit_p",
"MCTrack-hit proton; y; v_{1}", v1_nbx, -yvmax, yvmax, v1_nby, -1., 1.);
668 new TH2F(
"v1_rap_hit_pbar",
"MCTrack-hit antiproton;y; v_{1}", v1_nbx, -yvmax, yvmax, v1_nby, -1., 1.);
669 fa_v1_rap_hit_d =
new TH2F(
"v1_rap_hit_d",
"MCTrack-hit deuteron;y; v_{1}", v1_nbx, -yvmax, yvmax, v1_nby, -1., 1.);
670 fa_v1_rap_hit_t =
new TH2F(
"v1_rap_hit_t",
"MCTrack-hit triton; y; v_{1}", v1_nbx, -yvmax, yvmax, v1_nby, -1., 1.);
671 fa_v1_rap_hit_h =
new TH2F(
"v1_rap_hit_h",
"MCTrack-hit 3he; y; v_{1}", v1_nbx, -yvmax, yvmax, v1_nby, -1., 1.);
672 fa_v1_rap_hit_a =
new TH2F(
"v1_rap_hit_a",
"MCTrack-hit alpha; y; v_{1}", v1_nbx, -yvmax, yvmax, v1_nby, -1., 1.);
674 new TH2F(
"v1_rap_hit_imf",
"MCTrack-hit imf; y; v_{1}", v1_nbx, -yvmax, yvmax, v1_nby, -1., 1.);
677 new TH2F(
"v2_rap_hit_pip",
"MCTrack-hit pi-plus; y; v_{2}", v1_nbx, -yvmax, yvmax, v1_nby, -1., 1.);
679 new TH2F(
"v2_rap_hit_pim",
"MCTrack-hit pi-minus;y; v_{2}", v1_nbx, -yvmax, yvmax, v1_nby, -1., 1.);
680 fa_v2_rap_hit_kp =
new TH2F(
"v2_rap_hit_kp",
"MCTrack-hit k-plus; y; v_{2}", v1_nbx, -yvmax, yvmax, v1_nby, -1., 1.);
681 fa_v2_rap_hit_km =
new TH2F(
"v2_rap_hit_km",
"MCTrack-hit k-minus; y; v_{2}", v1_nbx, -yvmax, yvmax, v1_nby, -1., 1.);
682 fa_v2_rap_hit_p =
new TH2F(
"v2_rap_hit_p",
"MCTrack-hit proton; y; v_{2}", v1_nbx, -yvmax, yvmax, v1_nby, -1., 1.);
684 new TH2F(
"v2_rap_hit_pbar",
"MCTrack-hit antiproton;y; v_{2}", v1_nbx, -yvmax, yvmax, v1_nby, -1., 1.);
685 fa_v2_rap_hit_d =
new TH2F(
"v2_rap_hit_d",
"MCTrack-hit deuteron;y; v_{2}", v1_nbx, -yvmax, yvmax, v1_nby, -1., 1.);
686 fa_v2_rap_hit_t =
new TH2F(
"v2_rap_hit_t",
"MCTrack-hit triton; y; v_{2}", v1_nbx, -yvmax, yvmax, v1_nby, -1., 1.);
687 fa_v2_rap_hit_h =
new TH2F(
"v2_rap_hit_h",
"MCTrack-hit 3he; y; v_{2}", v1_nbx, -yvmax, yvmax, v1_nby, -1., 1.);
688 fa_v2_rap_hit_a =
new TH2F(
"v2_rap_hit_a",
"MCTrack-hit alpha; y; v_{2}", v1_nbx, -yvmax, yvmax, v1_nby, -1., 1.);
690 new TH2F(
"v2_rap_hit_imf",
"MCTrack-hit imf; y; v_{2}", v1_nbx, -yvmax, yvmax, v1_nby, -1., 1.);
694 fa_plab_sts_pip =
new TH1F(
"plab_sts_pip",
"MCTrack-sts pi-plus; p_{Lab}(GeV/c)", 100, 0., 10.);
695 fa_plab_sts_pim =
new TH1F(
"plab_sts_pim",
"MCTrack-sts pi-minus; p_{Lab}(GeV/c)", 100, 0., 10.);
696 fa_plab_sts_kp =
new TH1F(
"plab_sts_kp",
"MCTrack-sts k-plus; p_{Lab}(GeV/c)", 100, 0., 10.);
697 fa_plab_sts_km =
new TH1F(
"plab_sts_km",
"MCTrack-sts k-minus; p_{Lab}(GeV/c)", 100, 0., 10.);
698 fa_plab_sts_p =
new TH1F(
"plab_sts_p",
"MCTrack-sts proton; p_{Lab}(GeV/c)", 100, 0., 10.);
699 fa_plab_sts_pbar =
new TH1F(
"plab_sts_pbar",
"MCTrack-sts pbar; p_{Lab}(GeV/c)", 100, 0., 10.);
702 new TH2F(
"ptm_rap_sts_pip",
"MCTrack-sts pi-plus; y; p_{T}/m", ptm_nbx, ymin, ymax, ptm_nby, 0., ptmmax);
704 new TH2F(
"ptm_rap_sts_pim",
"MCTrack-sts pi-minus;y; p_{T}/m", ptm_nbx, ymin, ymax, ptm_nby, 0., ptmmax);
706 new TH2F(
"ptm_rap_sts_kp",
"MCTrack-sts k-plus; y; p_{T}/m", ptm_nbx, ymin, ymax, ptm_nby, 0., ptmmax);
708 new TH2F(
"ptm_rap_sts_km",
"MCTrack-sts k-minus; y; p_{T}/m", ptm_nbx, ymin, ymax, ptm_nby, 0., ptmmax);
710 new TH2F(
"ptm_rap_sts_p",
"MCTrack-sts proton; y; p_{T}/m", ptm_nbx, ymin, ymax, ptm_nby, 0., ptmmax);
712 new TH2F(
"ptm_rap_sts_pbar",
"MCTrack-sts antiproton;y;p_{T}/m", ptm_nbx, ymin, ymax, ptm_nby, 0., ptmmax);
714 new TH2F(
"ptm_rap_sts_d",
"MCTrack-sts deuteron;y; p_{T}/m", ptm_nbx, ymin, ymax, ptm_nby, 0., ptmmax);
716 new TH2F(
"ptm_rap_sts_t",
"MCTrack-sts triton; y; p_{T}/m", ptm_nbx, ymin, ymax, ptm_nby, 0., ptmmax);
718 new TH2F(
"ptm_rap_sts_h",
"MCTrack-sts 3he; y; p_{T}/m", ptm_nbx, ymin, ymax, ptm_nby, 0., ptmmax);
720 new TH2F(
"ptm_rap_sts_a",
"MCTrack-sts alpha; y; p_{T}/m", ptm_nbx, ymin, ymax, ptm_nby, 0., ptmmax);
722 new TH2F(
"ptm_rap_sts_imf",
"MCTrack-sts imf; y; p_{T}/m", ptm_nbx, ymin, ymax, ptm_nby, 0., ptmmax);
725 new TH2F(
"ptm_rap_glo_pip",
"MCTrack-glo pi-plus; y; p_{T}/m", ptm_nbx, ymin, ymax, ptm_nby, 0., ptmmax);
727 new TH2F(
"ptm_rap_glo_pim",
"MCTrack-glo pi-minus;y; p_{T}/m", ptm_nbx, ymin, ymax, ptm_nby, 0., ptmmax);
729 new TH2F(
"ptm_rap_glo_kp",
"MCTrack-glo k-plus; y; p_{T}/m", ptm_nbx, ymin, ymax, ptm_nby, 0., ptmmax);
731 new TH2F(
"ptm_rap_glo_km",
"MCTrack-glo k-minus; y; p_{T}/m", ptm_nbx, ymin, ymax, ptm_nby, 0., ptmmax);
733 new TH2F(
"ptm_rap_glo_p",
"MCTrack-glo proton; y; p_{T}/m", ptm_nbx, ymin, ymax, ptm_nby, 0., ptmmax);
735 new TH2F(
"ptm_rap_glo_pbar",
"MCTrack-glo antiproton;y;p_{T}/m", ptm_nbx, ymin, ymax, ptm_nby, 0., ptmmax);
737 new TH2F(
"ptm_rap_glo_d",
"MCTrack-glo deuteron;y; p_{T}/m", ptm_nbx, ymin, ymax, ptm_nby, 0., ptmmax);
739 new TH2F(
"ptm_rap_glo_t",
"MCTrack-glo triton; y; p_{T}/m", ptm_nbx, ymin, ymax, ptm_nby, 0., ptmmax);
741 new TH2F(
"ptm_rap_glo_h",
"MCTrack-glo 3he; y; p_{T}/m", ptm_nbx, ymin, ymax, ptm_nby, 0., ptmmax);
743 new TH2F(
"ptm_rap_glo_a",
"MCTrack-glo alpha; y; p_{T}/m", ptm_nbx, ymin, ymax, ptm_nby, 0., ptmmax);
745 new TH2F(
"ptm_rap_glo_imf",
"MCTrack-glo imf; y; p_{T}/m", ptm_nbx, ymin, ymax, ptm_nby, 0., ptmmax);
748 new TH2F(
"v1_rap_glo_pip",
"MCTrack-glo pi-plus; y; v_{1}", v1_nbx, -yvmax, yvmax, v1_nby, -1., 1.);
750 new TH2F(
"v1_rap_glo_pim",
"MCTrack-glo pi-minus;y; v_{1}", v1_nbx, -yvmax, yvmax, v1_nby, -1., 1.);
751 fa_v1_rap_glo_kp =
new TH2F(
"v1_rap_glo_kp",
"MCTrack-glo k-plus; y; v_{1}", v1_nbx, -yvmax, yvmax, v1_nby, -1., 1.);
752 fa_v1_rap_glo_km =
new TH2F(
"v1_rap_glo_km",
"MCTrack-glo k-minus; y; v_{1}", v1_nbx, -yvmax, yvmax, v1_nby, -1., 1.);
753 fa_v1_rap_glo_p =
new TH2F(
"v1_rap_glo_p",
"MCTrack-glo proton; y; v_{1}", v1_nbx, -yvmax, yvmax, v1_nby, -1., 1.);
755 new TH2F(
"v1_rap_glo_pbar",
"MCTrack-glo antiproton;y; v_{1}", v1_nbx, -yvmax, yvmax, v1_nby, -1., 1.);
756 fa_v1_rap_glo_d =
new TH2F(
"v1_rap_glo_d",
"MCTrack-glo deuteron;y; v_{1}", v1_nbx, -yvmax, yvmax, v1_nby, -1., 1.);
757 fa_v1_rap_glo_t =
new TH2F(
"v1_rap_glo_t",
"MCTrack-glo triton; y; v_{1}", v1_nbx, -yvmax, yvmax, v1_nby, -1., 1.);
758 fa_v1_rap_glo_h =
new TH2F(
"v1_rap_glo_h",
"MCTrack-glo 3he; y; v_{1}", v1_nbx, -yvmax, yvmax, v1_nby, -1., 1.);
759 fa_v1_rap_glo_a =
new TH2F(
"v1_rap_glo_a",
"MCTrack-glo alpha; y; v_{1}", v1_nbx, -yvmax, yvmax, v1_nby, -1., 1.);
761 new TH2F(
"v1_rap_glo_imf",
"MCTrack-glo imf; y; v_{1}", v1_nbx, -yvmax, yvmax, v1_nby, -1., 1.);
764 new TH2F(
"v2_rap_glo_pip",
"MCTrack-glo pi-plus; y; v_{2}", v1_nbx, -yvmax, yvmax, v1_nby, -1., 1.);
766 new TH2F(
"v2_rap_glo_pim",
"MCTrack-glo pi-minus;y; v_{2}", v1_nbx, -yvmax, yvmax, v1_nby, -1., 1.);
767 fa_v2_rap_glo_kp =
new TH2F(
"v2_rap_glo_kp",
"MCTrack-glo k-plus; y; v_{2}", v1_nbx, -yvmax, yvmax, v1_nby, -1., 1.);
768 fa_v2_rap_glo_km =
new TH2F(
"v2_rap_glo_km",
"MCTrack-glo k-minus; y; v_{2}", v1_nbx, -yvmax, yvmax, v1_nby, -1., 1.);
769 fa_v2_rap_glo_p =
new TH2F(
"v2_rap_glo_p",
"MCTrack-glo proton; y; v_{2}", v1_nbx, -yvmax, yvmax, v1_nby, -1., 1.);
771 new TH2F(
"v2_rap_glo_pbar",
"MCTrack-glo antiproton;y; v_{2}", v1_nbx, -yvmax, yvmax, v1_nby, -1., 1.);
772 fa_v2_rap_glo_d =
new TH2F(
"v2_rap_glo_d",
"MCTrack-glo deuteron;y; v_{2}", v1_nbx, -yvmax, yvmax, v1_nby, -1., 1.);
773 fa_v2_rap_glo_t =
new TH2F(
"v2_rap_glo_t",
"MCTrack-glo triton; y; v_{2}", v1_nbx, -yvmax, yvmax, v1_nby, -1., 1.);
774 fa_v2_rap_glo_h =
new TH2F(
"v2_rap_glo_h",
"MCTrack-glo 3he; y; v_{2}", v1_nbx, -yvmax, yvmax, v1_nby, -1., 1.);
775 fa_v2_rap_glo_a =
new TH2F(
"v2_rap_glo_a",
"MCTrack-glo alpha; y; v_{2}", v1_nbx, -yvmax, yvmax, v1_nby, -1., 1.);
777 new TH2F(
"v2_rap_glo_imf",
"MCTrack-glo imf; y; v_{2}", v1_nbx, -yvmax, yvmax, v1_nby, -1., 1.);
782 Float_t xrange = 750.;
783 Float_t yrange = 500.;
784 fwxy2 = nbinx * nbiny / 4. / xrange / yrange;
786 fa_xy_poi1 =
new TH2F(
"xy_poi1",
"TofPoint; ;", nbinx, -xrange, xrange, nbiny, -yrange, yrange);
787 fa_xy_poi2 =
new TH2F(
"xy_poi2",
"TofPoint /cm^{2}; ;", nbinx, -xrange, xrange, nbiny, -yrange, yrange);
788 fa_xy_poi3 =
new TH2F(
"xy_poi3",
"TofPoint /cm^{2}/s; ;", nbinx, -xrange, xrange, nbiny, -yrange, yrange);
790 fa_pv_poi =
new TH2F(
"pv_poi",
"TofPoint(all); velocity;momentum;", 100, 0., 32., 100., 0., 5.);
792 new TH2F(
"tm_poi",
"Tofpoi(all); momentum (GeV/c); M_{ToF}*sign(Z) (GeV/c^{2});", 100, 0., 10., 200., -1.5, 2.5);
793 fa_tm_poiprim =
new TH2F(
"tm_poiprim",
"Tofpoi(primary); momentum; Tofmass", 100, 0., 10., 200., -1.5, 2.5);
795 fa_xy_hit1 =
new TH2F(
"xy_hit1",
"TofHit; ;", nbinx, -xrange, xrange, nbiny, -yrange, yrange);
796 fa_xy_hit2 =
new TH2F(
"xy_hit2",
"TofHit /cm^{2}; ;", nbinx, -xrange, xrange, nbiny, -yrange, yrange);
797 fa_xy_hit3 =
new TH2F(
"xy_hit3",
"TofHit /cm^{2}/s; ;", nbinx, -xrange, xrange, nbiny, -yrange, yrange);
800 fa_tof_hit =
new TH1F(
"tof_hit",
"TofHit(all); t (ns); counts;", 100, 0., 100.);
801 fa_dtof_hit =
new TH1F(
"dtof_hit",
"TofHit(all); #Deltat (ns); counts;", 100, -1., 1.);
802 fa_tof_hitprim =
new TH1F(
"tof_hitprim",
"TofHit(prim); t (ns); counts;", 100, 0., 100.);
803 fa_pv_hit =
new TH2F(
"pv_hit",
"TofHit(all); velocity; momentum;", 100, 0., 32., 100., 0., 5.);
805 new TH2F(
"tm_hit",
"TofHit(all); momentum (GeV/c); M_{ToF}*sign(Z) (GeV/c^{2});", 100, 0., 10., 200., -1.5, 2.5);
806 fa_tm_hitprim =
new TH2F(
"tm_hitprim",
"TofHit(primary); momentum; Tofmass", 100, 0., 10., 200., -1.5, 2.5);
807 fa_tn_hit =
new TH1F(
"tn_hit",
"fastest TofHit(all); t (ns); counts;", 100, -1.0, 1.);
808 fa_t0_hit =
new TH1F(
"t0_hit",
"time zero; t0 (ns) ; counts;", 100, -0.5, 0.5);
809 fa_t0m_hit =
new TH1F(
"t0m_hit",
"average time zero; t0 (ns) ; counts;", 100, -0.1, 0.1);
810 fa_t0mn_hit =
new TH1F(
"t0mn_hit",
"average time zero hits; number of hits ; counts;", 100, 0., 500.);
811 fa_t0m_b_hit =
new TH2F(
"t0m_b_hit",
"average time zero; b (fm); t0 (ns) ; counts;", 28, 0., 14., 500, -0.3, 0.7);
813 new TH2F(
"t0mn_b_hit",
"average time zero hits; b(fm); number of hits ; counts;", 28, 0., 14., 100, 0., 500.);
815 fa_t0m_f_hit =
new TH1F(
"t0m_f_hit",
"average time zero forward; t0 (ns) ; counts;", 100, -0.1, 0.1);
816 fa_t0mn_f_hit =
new TH1F(
"t0mn_f_hit",
"average time zero hits forward; number of hits ; counts;", 100, 0., 500.);
818 new TH2F(
"t0m_f_b_hit",
"average time zero forward; b (fm); t0 (ns) ; counts;", 28, 0., 14., 500, -0.3, 0.7);
819 fa_t0mn_f_b_hit =
new TH2F(
"t0mn_f_b_hit",
"average time zero hits forward ; b(fm); number of hits ; counts;", 28, 0.,
822 fa_t0m_nf_hit =
new TH1F(
"t0m_nf_hit",
"average time zero noforward; t0 (ns) ; counts;", 100, -0.1, 0.1);
823 fa_t0mn_nf_hit =
new TH1F(
"t0mn_nf_hit",
"average time zero hits noforward; number of hits ; counts;", 100, 0., 500.);
825 new TH2F(
"t0m_nf_b_hit",
"average time zero noforward; b (fm); t0 (ns) ; counts;", 28, 0., 14., 500, -0.3, 0.7);
826 fa_t0mn_nf_b_hit =
new TH2F(
"t0mn_nf_b_hit",
"average time zero hits noforward ; b(fm); number of hits ; counts;", 28,
827 0., 14., 100, 0., 500.);
829 fa_dxx =
new TH2F(
"dxx",
"TofHit - TofPoint; x (cm); #Delta x (cm);", 100, -xrange, xrange, 50., -10., 10.);
830 fa_dxy =
new TH2F(
"dxy",
"TofHit - TofPoint; y (cm); #Delta x (cm);", 100, -yrange, yrange, 50., -10., 10.);
831 fa_dxz =
new TH2F(
"dxz",
"TofHit - TofPoint; z (cm); #Delta x (cm);", 100, 400., 650., 50., -10., 10.);
832 fa_dyx =
new TH2F(
"dyx",
"TofHit - TofPoint; x (cm); #Delta y (cm);", 100, -xrange, xrange, 50., -10., 10.);
833 fa_dyy =
new TH2F(
"dyy",
"TofHit - TofPoint; y (cm); #Delta y (cm);", 100, -yrange, yrange, 50., -10., 10.);
834 fa_dyz =
new TH2F(
"dyz",
"TofHit - TofPoint; z (cm); #Delta y (cm);", 100, 400., 650., 50., -10., 10.);
835 fa_dzx =
new TH2F(
"dzx",
"TofHit - TofPoint; x (cm); #Delta z (cm);", 100, -xrange, xrange, 50, -20., 20.);
836 fa_dzy =
new TH2F(
"dzy",
"TofHit - TofPoint; y (cm); #Delta z (cm);", 100, -yrange, yrange, 50, -20., 20.);
837 fa_dzz =
new TH2F(
"dzz",
"TofHit - TofPoint; z (cm); #Delta z (cm);", 100, 400., 650., 50, -20., 20.);
839 fa_hit_ch =
new TH1F(
"hit_ch",
"TofHits; channel; rate (Hz/s);", 50000, 0., 50000.);
840 fa_dhit_ch =
new TH2F(
"dhit_ch",
"Tof Double Hits; channel; counts;", 50000, 0., 50000., 2, 0.5, 2.5);
842 fa_xy_glo1 =
new TH2F(
"xy_glo1",
"GlobalTrack(all); x (cm); y(cm);", nbinx, -xrange, xrange, nbiny, -yrange, yrange);
844 new TH2F(
"xy_glo_pip",
"GlobalTrack(pip); x (cm); y(cm);", nbinx, -xrange, xrange, nbiny, -yrange, yrange);
846 new TH2F(
"xy_glo_pim",
"GlobalTrack(pim); x (cm); y(cm);", nbinx, -xrange, xrange, nbiny, -yrange, yrange);
847 fa_xy_glo_p =
new TH2F(
"xy_glo_p",
"GlobalTrack(p); x (cm); y(cm);", nbinx, -xrange, xrange, nbiny, -yrange, yrange);
849 new TH2F(
"xy_glo_pbar",
"GlobalTrack(pbar); x (cm); y(cm);", nbinx, -xrange, xrange, nbiny, -yrange, yrange);
851 new TH2F(
"xy_glo_kp",
"GlobalTrack(kp); x (cm); y(cm);", nbinx, -xrange, xrange, nbiny, -yrange, yrange);
853 new TH2F(
"xy_glo_km",
"GlobalTrack(km); x (cm); y(cm);", nbinx, -xrange, xrange, nbiny, -yrange, yrange);
854 fa_xy_glo_d =
new TH2F(
"xy_glo_d",
"GlobalTrack(d); x (cm); y(cm);", nbinx, -xrange, xrange, nbiny, -yrange, yrange);
855 fa_xy_glo_t =
new TH2F(
"xy_glo_t",
"GlobalTrack(t); x (cm); y(cm);", nbinx, -xrange, xrange, nbiny, -yrange, yrange);
856 fa_xy_glo_h =
new TH2F(
"xy_glo_h",
"GlobalTrack(h); x (cm); y(cm);", nbinx, -xrange, xrange, nbiny, -yrange, yrange);
857 fa_xy_glo_a =
new TH2F(
"xy_glo_a",
"GlobalTrack(a); x (cm); y(cm);", nbinx, -xrange, xrange, nbiny, -yrange, yrange);
859 fa_TofTrackMul =
new TH1F(
"TofTrackMul",
"number of assigned TofTrack / global track; multiplicity; ", 50, 0., 50.);
860 fa_VtxB =
new TH1F(
"VtxB",
"Chi2 to primary vertex; ", 100, 0., 20.);
862 Float_t TMMIN = -1.5;
865 fa_pv_glo =
new TH2F(
"pv_glo",
"GlobalTrack(all); velocity; momentum;", 100, 0., 32., 100., 0., 5.);
866 fa_tm_glo =
new TH2F(
"tm_glo",
"GlobalTrack(all); momentum (GeV/c); M_{ToF}*sign(Z) (GeV/c^{2});", 100, 0., 10.,
867 TMYBIN, TMMIN, TMMAX);
868 fa_chi2_mom_glo =
new TH2F(
"chi2_mom_glo",
"GlobalTrack(all); momentum; chi2;", 100, 0., 10., 100., 0., 50.);
870 new TH2F(
"chi2_mom_gloprim",
"GlobalTrack(primaries); momentum; chi2;", 100, 0., 10., 100., 0., 50.);
871 fa_len_mom_glo =
new TH2F(
"len_mom_glo",
"GlobalTrack(all); momentum; len;", 100, 0., 10., 100., 0., 1500.);
872 fa_tm_gloprim =
new TH2F(
"tm_gloprim",
"GlobalTrack(primary); momentum; Tofmass", 100, 0., 10., TMYBIN, TMMIN, TMMAX);
873 fa_tm_glovtxb =
new TH2F(
"tm_glovtxb",
"GlobalTrack(vtxb); momentum (GeV/c); M_{ToF}*sign(Z) (GeV/c^{2});", 100, 0.,
874 10., TMYBIN, TMMIN, TMMAX);
876 new TH2F(
"tm_gloprimvtxb",
"GlobalTrack(prim,vtxb); momentum; Tofmass", 100, 0., 10., TMYBIN, TMMIN, TMMAX);
877 fa_tm_glomis =
new TH2F(
"tm_glomis",
"GlobalTrack(mis); momentum (GeV/c); M_{ToF}*sign(Z) (GeV/c^{2});", 100, 0., 10.,
878 TMYBIN, TMMIN, TMMAX);
879 fa_tm_glo_pip =
new TH2F(
"tm_glo_pip",
"GlobalTrack(pip); momentum (GeV/c); M_{ToF}*sign(Z) (GeV/c^{2});", 100, 0.,
880 10., TMYBIN, TMMIN, TMMAX);
881 fa_tm_glo_pim =
new TH2F(
"tm_glo_pim",
"GlobalTrack(pim); momentum (GeV/c); M_{ToF}*sign(Z) (GeV/c^{2});", 100, 0.,
882 10., TMYBIN, TMMIN, TMMAX);
883 fa_tm_glo_kp =
new TH2F(
"tm_glo_kp",
"GlobalTrack(kp); momentum (GeV/c); M_{ToF}*sign(Z) (GeV/c^{2});", 100, 0., 10.,
884 TMYBIN, TMMIN, TMMAX);
885 fa_tm_glo_km =
new TH2F(
"tm_glo_km",
"GlobalTrack(km); momentum (GeV/c); M_{ToF}*sign(Z) (GeV/c^{2});", 100, 0., 10.,
886 TMYBIN, TMMIN, TMMAX);
887 fa_tm_glo_p =
new TH2F(
"tm_glo_p",
"GlobalTrack(p); momentum (GeV/c); M_{ToF}*sign(Z) (GeV/c^{2});", 100, 0., 10.,
888 TMYBIN, TMMIN, TMMAX);
889 fa_tm_glo_pbar =
new TH2F(
"tm_glo_pbar",
"GlobalTrack(pbar); momentum (GeV/c); M_{ToF}*sign(Z) (GeV/c^{2});", 100, 0.,
890 10., TMYBIN, TMMIN, TMMAX);
891 fa_tm_glo_d =
new TH2F(
"tm_glo_d",
"GlobalTrack(d); momentum (GeV/c); M_{ToF}*sign(Z) (GeV/c^{2});", 100, 0., 10.,
892 TMYBIN, TMMIN, TMMAX);
893 fa_tm_glo_t =
new TH2F(
"tm_glo_t",
"GlobalTrack(t); momentum (GeV/c); M_{ToF}*sign(Z) (GeV/c^{2});", 100, 0., 10.,
894 TMYBIN, TMMIN, TMMAX);
895 fa_tm_glo_h =
new TH2F(
"tm_glo_h",
"GlobalTrack(h); momentum (GeV/c); M_{ToF}*sign(Z) (GeV/c^{2});", 100, 0., 10.,
896 TMYBIN, TMMIN, TMMAX);
897 fa_tm_glo_a =
new TH2F(
"tm_glo_a",
"GlobalTrack(a); momentum (GeV/c); M_{ToF}*sign(Z) (GeV/c^{2});", 100, 0., 10.,
898 TMYBIN, TMMIN, TMMAX);
900 Double_t M2MIN = -0.4;
901 Double_t M2MAX = 1.4;
903 fa_m2mom_glo =
new TH2F(
"m2mom_glo",
"GlobalTrack(all); p ^{.} sign(Z)(GeV); M_{ToF}^{2} (GeV^{2});", 200, -10., 10.,
904 M2YBIN, M2MIN, M2MAX);
905 fa_m2mom_glovtxb =
new TH2F(
"m2mom_glovtxb",
"GlobalTrack(vtxb); p ^{.} sign(Z)(GeV); M_{ToF}^{2} (GeV^{2});", 200,
906 -10., 10., M2YBIN, M2MIN, M2MAX);
907 fa_m2mom_gloprim =
new TH2F(
"m2mom_gloprim",
"GlobalTrack(prim); p ^{.} sign(Z)(GeV); M_{ToF}^{2} (GeV^{2});", 200,
908 -10., 10., M2YBIN, M2MIN, M2MAX);
910 new TH2F(
"m2mom_gloprimvtxb",
"GlobalTrack(primvtxb); p ^{.} sign(Z)(GeV); M_{ToF}^{2} (GeV^{2});", 200, -10., 10.,
911 M2YBIN, M2MIN, M2MAX);
912 fa_m2mom_glo_pip =
new TH2F(
"m2mom_glo_pip",
"GlobalTrack(pip); p ^{.} sign(Z)(GeV); M_{ToF}^{2} (GeV^{2});", 200,
913 -10., 10., M2YBIN, M2MIN, M2MAX);
914 fa_m2mom_glo_pim =
new TH2F(
"m2mom_glo_pim",
"GlobalTrack(pim); p ^{.} sign(Z)(GeV); M_{ToF}^{2} (GeV^{2});", 200,
915 -10., 10., M2YBIN, M2MIN, M2MAX);
916 fa_m2mom_glo_kp =
new TH2F(
"m2mom_glo_kp",
"GlobalTrack(kp); p ^{.} sign(Z)(GeV); M_{ToF}^{2} (GeV^{2});", 200, -10.,
917 10., M2YBIN, M2MIN, M2MAX);
918 fa_m2mom_glo_km =
new TH2F(
"m2mom_glo_km",
"GlobalTrack(km); p ^{.} sign(Z)(GeV); M_{ToF}^{2} (GeV^{2});", 200, -10.,
919 10., M2YBIN, M2MIN, M2MAX);
920 fa_m2mom_glo_p =
new TH2F(
"m2mom_glo_p",
"GlobalTrack(p); p ^{.} sign(Z)(GeV); M_{ToF}^{2} (GeV^{2});", 200, -10.,
921 10., M2YBIN, M2MIN, M2MAX);
922 fa_m2mom_glo_pbar =
new TH2F(
"m2mom_glo_pbar",
"GlobalTrack(pbar); p ^{.} sign(Z)(GeV); M_{ToF}^{2} (GeV^{2});", 200,
923 -10., 10., M2YBIN, M2MIN, M2MAX);
925 fa_pMCmom_glo =
new TH2F(
"pMCmom_glo",
"GlobalTrack(all); momentum; p_{MC};", 100, 0., 10., 100, 0., 10.);
926 fa_pMCmom_glo_pip =
new TH2F(
"pMCmom_glo_pip",
"GlobalTrack(pip); momentum; p_{MC};", 100, 0., 10., 100, 0., 10.);
927 fa_pMCmom_glo_pim =
new TH2F(
"pMCmom_glo_pim",
"GlobalTrack(pim); momentum; p_{MC};", 100, 0., 10., 100, 0., 10.);
928 fa_pMCmom_glo_kp =
new TH2F(
"pMCmom_glo_kp",
"GlobalTrack(kp); momentum; p_{MC};", 100, 0., 10., 100, 0., 10.);
929 fa_pMCmom_glo_km =
new TH2F(
"pMCmom_glo_km",
"GlobalTrack(km); momentum; p_{MC};", 100, 0., 10., 100, 0., 10.);
930 fa_pMCmom_glo_p =
new TH2F(
"pMCmom_glo_p",
"GlobalTrack(p); momentum; p_{MC};", 100, 0., 10., 100, 0., 10.);
931 fa_pMCmom_glo_pbar =
new TH2F(
"pMCmom_glo_pbar",
"GlobalTrack(pbar); momentum; p_{MC};", 100, 0., 10., 100, 0., 10.);
932 fa_pMCmom_glo_d =
new TH2F(
"pMCmom_glo_d",
"GlobalTrack(d); momentum; p_{MC};", 100, 0., 10., 100, 0., 10.);
933 fa_pMCmom_glo_t =
new TH2F(
"pMCmom_glo_t",
"GlobalTrack(t); momentum; p_{MC};", 100, 0., 10., 100, 0., 10.);
934 fa_pMCmom_glo_h =
new TH2F(
"pMCmom_glo_h",
"GlobalTrack(h); momentum; p_{MC};", 100, 0., 10., 100, 0., 10.);
935 fa_pMCmom_glo_a =
new TH2F(
"pMCmom_glo_a",
"GlobalTrack(a); momentum; p_{MC};", 100, 0., 10., 100, 0., 10.);
937 Double_t LenDisMax = 20.;
939 new TH2F(
"LenDismom_glo",
"GlobalTrack(all); momentum; len-dis;", 100, 0., 10., 100., -LenDisMax, LenDisMax);
941 new TH2F(
"LenDismom_glo_pip",
"GlobalTrack(pip); momentum; len-dis;", 100, 0., 10., 100., -LenDisMax, LenDisMax);
943 new TH2F(
"LenDismom_glo_pim",
"GlobalTrack(pim); momentum; len-dis;", 100, 0., 10., 100., -LenDisMax, LenDisMax);
945 new TH2F(
"LenDismom_glo_kp",
"GlobalTrack(kp); momentum; len-dis;", 100, 0., 10., 100., -LenDisMax, LenDisMax);
947 new TH2F(
"LenDismom_glo_km",
"GlobalTrack(km); momentum; len-dis;", 100, 0., 10., 100., -LenDisMax, LenDisMax);
949 new TH2F(
"LenDismom_glo_p",
"GlobalTrack(p); momentum; len-dis;", 100, 0., 10., 100., -LenDisMax, LenDisMax);
951 new TH2F(
"LenDismom_glo_pbar",
"GlobalTrack(pbar); momentum; len-dis;", 100, 0., 10., 100., -LenDisMax, LenDisMax);
953 new TH2F(
"LenDismom_glo_d",
"GlobalTrack(d); momentum; len-dis;", 100, 0., 10., 100., -LenDisMax, LenDisMax);
955 new TH2F(
"LenDismom_glo_t",
"GlobalTrack(t); momentum; len-dis;", 100, 0., 10., 100., -LenDisMax, LenDisMax);
957 new TH2F(
"LenDismom_glo_h",
"GlobalTrack(h); momentum; len-dis;", 100, 0., 10., 100., -LenDisMax, LenDisMax);
959 new TH2F(
"LenDismom_glo_a",
"GlobalTrack(a); momentum; len-dis;", 100, 0., 10., 100., -LenDisMax, LenDisMax);
961 fa_LenMcLenGlomom_glo =
new TH2F(
"LenMcLenGlomom_glo",
"GlobalTrack(all); momentum [GeV]; len glo - len MC pnt[cm];",
962 100, 0., 10., 400., -LenDisMax, LenDisMax);
963 fa_LenMcDismom_glo =
new TH2F(
"LenMcDismom_glo",
"GlobalTrack(all); momentum [GeV]; len MC pnt - dis [cm];", 100, 0.,
964 10., 400., -LenDisMax, LenDisMax);
968 fa_w_mom_glo =
new TH2F(
"w_mom_glo",
"GlobalTrack(all); momentum; weight;", 100, 0., 10., WYBIN, 0., WYMAX);
969 fa_w_mom_glo_pip =
new TH2F(
"w_mom_glo_pip",
"GlobalTrack(pip); momentum; weight;", 100, 0., 10., WYBIN, 0., WYMAX);
970 fa_w_mom_glo_pim =
new TH2F(
"w_mom_glo_pim",
"GlobalTrack(pim); momentum; weight;", 100, 0., 10., WYBIN, 0., WYMAX);
971 fa_w_mom_glo_kp =
new TH2F(
"w_mom_glo_kp",
"GlobalTrack(kp); momentum; weight;", 100, 0., 10., WYBIN, 0., WYMAX);
972 fa_w_mom_glo_km =
new TH2F(
"w_mom_glo_km",
"GlobalTrack(km); momentum; weight;", 100, 0., 10., WYBIN, 0., WYMAX);
973 fa_w_mom_glo_p =
new TH2F(
"w_mom_glo_p",
"GlobalTrack(p); momentum; weight;", 100, 0., 10., WYBIN, 0., WYMAX);
975 new TH2F(
"w_mom_glo_pbar",
"GlobalTrack(pbar); momentum; weight;", 100, 0., 10., WYBIN, 0., WYMAX);
976 fa_w_mom_glo_d =
new TH2F(
"w_mom_glo_d",
"GlobalTrack(d); momentum; weight;", 100, 0., 10., WYBIN, 0., WYMAX);
977 fa_w_mom_glo_t =
new TH2F(
"w_mom_glo_t",
"GlobalTrack(t); momentum; weight;", 100, 0., 10., WYBIN, 0., WYMAX);
978 fa_w_mom_glo_h =
new TH2F(
"w_mom_glo_h",
"GlobalTrack(h); momentum; weight;", 100, 0., 10., WYBIN, 0., WYMAX);
979 fa_w_mom_glo_a =
new TH2F(
"w_mom_glo_a",
"GlobalTrack(a); momentum; weight;", 100, 0., 10., WYBIN, 0., WYMAX);
980 fa_w_mom_glomis =
new TH2F(
"w_mom_glomis",
"GlobalTrack(mis); momentum; weight;", 100, 0., 10., WYBIN, 0., WYMAX);
984 new TH2F(
"mul_b_gen",
"Centrality - gen;impact parameter b; multiplicity;", 15, 0., 15., 100, 0., 10000.);
986 new TH2F(
"mul_b_poi",
"Centrality - poi;impact parameter b; multiplicity;", 15, 0., 15., 100, 0., 10000.);
988 new TH2F(
"mul_b_hit",
"Centrality - hit;impact parameter b; multiplicity;", 15, 0., 15., 100, 0., 2000.);
990 new TH2F(
"mul_b_glo",
"Centrality - glo;impact parameter b; multiplicity;", 15, 0., 15., 100, 0., 2000.);
992 new TH2F(
"mul_b_had",
"Centrality - had;impact parameter b; multiplicity;", 15, 0., 15., 100, 0., 2000.);
995 fa_cdrp_b_gen =
new TH2F(
"cdrp_b_gen",
"reaction plane resolution - gen;impact parameter b; cos#Delta#phi_{rp};", 15,
996 0., 15., 20, -1., 1.);
998 new TH2F(
"drp_b_gen",
"#Delta#phi-reaction plane - gen ;impact parameter b;", 15, 0., 15., 36, -180., 180.);
1000 new TH2F(
"phirp_b_gen",
"#phi_{reaction plane} - gen ;impact parameter b;", 15, 0., 15., 36, -180., 180.);
1002 new TH2F(
"phgrp_b_gen",
"#phi_{G reaction plane} - gen ;impact parameter b;", 15, 0., 15., 36, -180., 180.);
1003 fa_phphrp_gen =
new TH2F(
"phphrp_gen",
"#phi_#phi - gen ;#phi_{rec}; #phi_{gen} ;", 36, -180., 180., 36, -180., 180.);
1005 "#Delta#phi_{G}-reaction plane - gen ;impact "
1006 "parameter b;#phi_{rec}-#phi_{gen}",
1007 15, 0., 15., 36, -180., 180.);
1009 "#Delta#phi_{G}-reaction plane - poi ;impact "
1010 "parameter b;#phi_{rec}-#phi_{gen}",
1011 15, 0., 15., 36, -180., 180.);
1013 "#Delta#phi_{G}-reaction plane - hit ;impact "
1014 "parameter b;#phi_{rec}-#phi_{gen}",
1015 15, 0., 15., 36, -180., 180.);
1017 "#Delta#phi_{G}-reaction plane - glo ;impact "
1018 "parameter b;#phi_{rec}-#phi_{gen}",
1019 15, 0., 15., 36, -180., 180.);
1022 "reaction plane resolution - gen;impact parameter "
1023 "b;cos(#phi_{rec}-#phi_{gen})",
1024 15, 0., 15., 20, -1., 1.);
1026 "reaction plane resolution - poi;impact parameter "
1027 "b;cos(#phi_{rec}-#phi_{gen})",
1028 15, 0., 15., 20, -1., 1.);
1030 "reaction plane resolution - hit;impact parameter "
1031 "b;cos(#phi_{rec}-#phi_{gen})",
1032 15, 0., 15., 20, -1., 1.);
1034 "reaction plane resolution - glo;impact parameter "
1035 "b;cos(#phi_{rec}-#phi_{gen})",
1036 15, 0., 15., 20, -1., 1.);
1038 "reaction plane resolution - had;impact parameter "
1039 "b;cos(#phi_{rec}-#phi_{gen})",
1040 15, 0., 15., 20, -1., 1.);
1042 fa_cdrp_b_poi =
new TH2F(
"cdrp_b_poi",
"reaction plane resolution - poi;impact parameter b; cos#Delta#phi_{rp};", 15,
1043 0., 15., 20, -1., 1.);
1045 new TH2F(
"drp_b_poi",
"#Delta#phi-reaction plane - poi ;impact parameter b;", 15, 0., 15., 36, -180., 180.);
1047 fa_cdrp_b_hit =
new TH2F(
"cdrp_b_hit",
"reaction plane resolution - hit;impact parameter b; cos#Delta#phi_{rp};", 15,
1048 0., 15., 20, -1., 1.);
1050 new TH2F(
"drp_b_hit",
"#Delta#phi-reaction plane - hit ;impact parameter b;", 15, 0., 15., 36, -180., 180.);
1052 fa_cdrp_b_glo =
new TH2F(
"cdrp_b_glo",
"reaction plane resolution - glo;impact parameter b; cos#Delta#phi_{rp};", 15,
1053 0., 15., 20, -1., 1.);
1055 new TH2F(
"drp_b_glo",
"#Delta#phi-reaction plane - glo ;impact parameter b;", 15, 0., 15., 36, -180., 180.);
1057 fa_cdrp_b_had =
new TH2F(
"cdrp_b_had",
"reaction plane resolution - had;impact parameter b; cos#Delta#phi_{rp};", 15,
1058 0., 15., 20, -1., 1.);
1060 new TH2F(
"drp_b_had",
"#Delta#phi-reaction plane - had ;impact parameter b;", 15, 0., 15., 36, -180., 180.);
1062 fa_phirp_gen =
new TH1F(
"phirp_gen",
"#phi_{reaction plane} - gen ;#phi_{RPgen};", 72, -180., 180.);
1063 fa_phirp_poi =
new TH1F(
"phirp_poi",
"#phi_{reaction plane} - poi ;#phi_{RPpoi};", 72, -180., 180.);
1064 fa_phirp_hit =
new TH1F(
"phirp_hit",
"#phi_{reaction plane} - hit ;#phi_{RPhit};", 72, -180., 180.);
1065 fa_phirp_glo =
new TH1F(
"phirp_glo",
"#phi_{reaction plane} - glo ;#phi_{RPglo};", 72, -180., 180.);
1066 fa_phirp_had =
new TH1F(
"phirp_had",
"#phi_{reaction plane} - had ;#phi_{RPhad};", 72, -180., 180.);
1068 fa_phirps_gen =
new TH1F(
"phirps_gen",
"#phi_{reaction plane} - gen ;#phi_{sRPgen};", 72, -180., 180.);
1069 fa_phirps_poi =
new TH1F(
"phirps_poi",
"#phi_{reaction plane} - poi ;#phi_{sRPpoi};", 72, -180., 180.);
1070 fa_phirps_hit =
new TH1F(
"phirps_hit",
"#phi_{reaction plane} - hit ;#phi_{sRPhit};", 72, -180., 180.);
1071 fa_phirps_glo =
new TH1F(
"phirps_glo",
"#phi_{reaction plane} - glo ;#phi_{sRPglo};", 72, -180., 180.);
1073 fhwdist =
new TH2F(
"hwdist",
"matching wdist; p (GeV/c); dist;", 100, 0., 10., 50, 0., 10.);
1074 fhwmindelmass =
new TH2F(
"hwmindelmass",
"matching wmindelmass ; p (GeV/c); #Delta M;", 100, 0., 10., 50, 0., 1.);
1075 fhwminlen =
new TH2F(
"hwminlen",
"matching wminlen ; p (GeV/c); MinPathlength;", 100, 0., 10., 50, 0., 2.);
1076 fhwdelp =
new TH2F(
"hwdelp",
"matching wdelp ; p (GeV/c); #Delta p/p;", 100, 0., 10., 50, 0., 1.);
1078 fhTofTrkDx =
new TH1F(
"hTofTrkDx",
" x - position resolution ; #deltax;", 50, -1., 1.);
1079 fhTofTrkDy =
new TH1F(
"hTofTrkDy",
" y - position resolution ; #deltay;", 50, -1., 1.);
1080 fhTofTrkDxsel =
new TH1F(
"hTofTrkDxsel",
" x - position resolution ; #deltax;", 50, -1., 1.);
1081 fhTofTrkDysel =
new TH1F(
"hTofTrkDysel",
" y - position resolution ; #deltay;", 50, -1., 1.);
1083 cout <<
"CbmHadronAnalysis::CreateHistogramms: histograms booked in directory " << gDirectory->GetName() << endl;
1376 Int_t pdgCode, Np1, Np2;
1377 Float_t Qx1, Qy1, Qx2, Qy2, phirp1, phirp2, phirp, delrp, rp_weight;
1378 Float_t RADDEG = 57.29577951;
1379 Float_t p_MC, px_MC, py_MC, pz_MC;
1381 Float_t MaxT0 = 0.1;
1383 Int_t TrackP[100000];
1385 Bool_t use_pions_for_flow = kTRUE;
1409 LOG(debug) <<
"<D> HadronAnalysis::Exec starting with MCtrks " <<
nMCTracks <<
", TofPoi " <<
nTofPoints
1415 if (fair::Logger::Logging(fair::Severity::debug)) {
1416 for (Int_t j = 0; j <
nTofHits; j++) {
1418 if (NULL == TofHit)
continue;
1419 LOG(debug) << Form(
"TofHit %d, addr 0x%08x, x %6.1f, y %6.1f, z %6.1f, t %6.1f ", j, TofHit->GetAddress(),
1420 TofHit->GetX(), TofHit->GetY(), TofHit->GetZ(), TofHit->GetTime());
1425 for (Int_t j = 0; j <
nTofHits; j++) {
1427 if (NULL == TofHit)
continue;
1428 if (TofHit->GetZ() == 0.) dT0 = TofHit->GetTime();
1432 for (Int_t j = 0; j <
nTofHits; j++) {
1434 if (NULL == TofHit)
continue;
1435 TofHit->SetTime(TofHit->GetTime() - dT0);
1439 if (fair::Logger::Logging(fair::Severity::debug)) {
1440 for (Int_t j = 0; j <
nTofHits; j++) {
1442 if (NULL == TofHit)
continue;
1443 LOG(debug) << Form(
"TofHit %d, addr 0x%08x, x %6.1f, y %6.1f, z %6.1f, t %6.1f ", j, TofHit->GetAddress(),
1444 TofHit->GetX(), TofHit->GetY(), TofHit->GetZ(), TofHit->GetTime());
1457 Int_t IndTofTrack_TofHit[
nTofHits][MAXNHT];
1478 if (pdgCode < 100000000) {
1479 if (211 != TMath::Abs(pdgCode) &&
1480 321 != TMath::Abs(pdgCode) &&
1481 2212 != TMath::Abs(pdgCode))
1483 LOG(info) <<
"skip pdgCode " << pdgCode;
1488 mfrag = (pdgCode % 1000) / 10 * .931494028;
1492 Float_t Phip = RADDEG * atan2(MCTrack->
GetPy(), MCTrack->
GetPx());
1494 if (dphi < -180.) { dphi += 360.; };
1495 if (dphi > 180.) { dphi -= 360.; };
1496 dphi = dphi / RADDEG;
1499 LOG(info) <<
" Track k=" << k <<
", pdgCode = " << pdgCode <<
" Mass " << MCTrack->
GetMass() <<
"," << mfrag
1508 if (use_pions_for_flow && TMath::Abs((MCTrack->
GetRapidity() - yrp_mid) / yrp_mid) >
GetDY()
1527 if (use_pions_for_flow && TMath::Abs((MCTrack->
GetRapidity() - yrp_mid) / yrp_mid) >
GetDY()
1666 if (rp_weight != 0.) {
1667 if (gRandom->Uniform(1) > 0.5) {
1669 Qx1 = Qx1 + rp_weight * MCTrack->
GetPx();
1670 Qy1 = Qy1 + rp_weight * MCTrack->
GetPy();
1674 Qx2 = Qx2 + rp_weight * MCTrack->
GetPx();
1675 Qy2 = Qy2 + rp_weight * MCTrack->
GetPy();
1680 if (Np1 > 0 && Np2 > 0) {
1681 phirp1 = atan2(Qy1, Qx1);
1682 phirp2 = atan2(Qy2, Qx2);
1684 TH1F* phirp_gen_fpar =
fflowFile->Get<TH1F>(
"phirps_gen_fpar");
1686 for (
int j = 0; j < 4; j++) {
1687 Float_t i = (float) (j + 1);
1688 dphir += (-phirp_gen_fpar->GetBinContent(j) * TMath::Cos(i * phirp1)
1689 + phirp_gen_fpar->GetBinContent(j + 4) * TMath::Sin(i * phirp1))
1695 for (
int j = 0; j < 4; j++) {
1696 Float_t i = (float) (j + 1);
1697 dphir += (-phirp_gen_fpar->GetBinContent(j) * TMath::Cos(i * phirp2)
1698 + phirp_gen_fpar->GetBinContent(j + 4) * TMath::Sin(i * phirp2))
1704 delrp = (phirp1 - phirp2);
1708 cout <<
"<D-gen> Impact parameter " <<
fMCEventHeader->GetB() <<
", delrp = " << delrp << endl;
1711 delrp = delrp * RADDEG;
1712 while (delrp < -180.) {
1715 while (delrp > 180.) {
1720 phirp = RADDEG * atan2(Qy1 + Qy2, Qx1 + Qx2);
1721 while (phirp < -180.) {
1724 while (phirp > 180.) {
1728 TH1F* phirp_gen_fpar =
fflowFile->Get<TH1F>(
"phirp_gen_fpar");
1730 for (
int j = 0; j < 4; j++) {
1731 Float_t i = (float) (j + 1);
1734 dphir += ((-phirp_gen_fpar->GetBinContent(j) * TMath::Cos(i * phirp / RADDEG)
1735 + phirp_gen_fpar->GetBinContent(j + 4) * TMath::Sin(i * phirp / RADDEG))
1740 phirp += dphir * RADDEG;
1741 while (phirp < -180.) {
1744 while (phirp > 180.) {
1749 while (delrp < -180.) {
1752 while (delrp > 180.) {
1779 if (NULL == TofPoint) {
1780 LOG(warning) <<
" Missing TofPoint at " << l <<
", mul " <<
nTofPoints;
1783 Int_t k = TofPoint->GetTrackID();
1786 if (pdgCode > 100000000) {
1787 mfrag = (pdgCode % 1000) / 10 * .931494028;
1789 px_MC = MCTrack->
GetPx();
1790 py_MC = MCTrack->
GetPy();
1791 pz_MC = MCTrack->
GetPz();
1792 p_MC =
sqrt(px_MC * px_MC + py_MC * py_MC + pz_MC * pz_MC);
1794 if (TrackP[k] == 0) {
1797 fa_xy_poi1->Fill(TofPoint->GetX(), TofPoint->GetY());
1800 Float_t vel = TofPoint->GetLength() / TofPoint->GetTime();
1801 Float_t bet = vel /
clight;
1802 if (bet > 0.99999) { bet = 0.99999; }
1803 Float_t tofmass = p_MC / bet * TMath::Sqrt(1. - bet * bet) * TMath::Sign(1, pdgCode);
1811 Float_t Phip = RADDEG * atan2(MCTrack->
GetPy(), MCTrack->
GetPx());
1813 if (dphi < -180.) { dphi += 360.; };
1814 if (dphi > 180.) { dphi -= 360.; };
1815 dphi = dphi / RADDEG;
1821 if (use_pions_for_flow && TMath::Abs((MCTrack->
GetRapidity() - yrp_mid) / yrp_mid) >
GetDY()
1837 if (use_pions_for_flow && TMath::Abs((MCTrack->
GetRapidity() - yrp_mid) / yrp_mid) >
GetDY()
1966 if (rp_weight != 0.) {
1967 if (gRandom->Uniform(1) > 0.5) {
1969 Qx1 = Qx1 + rp_weight * MCTrack->
GetPx();
1970 Qy1 = Qy1 + rp_weight * MCTrack->
GetPy();
1974 Qx2 = Qx2 + rp_weight * MCTrack->
GetPx();
1975 Qy2 = Qy2 + rp_weight * MCTrack->
GetPy();
1980 if (Np1 > 0 && Np2 > 0) {
1981 phirp1 = atan2(Qy1, Qx1);
1982 phirp2 = atan2(Qy2, Qx2);
1984 TH1F* phirp_poi_fpar =
fflowFile->Get<TH1F>(
"phirps_poi_fpar");
1986 for (
int j = 0; j < 4; j++) {
1987 Float_t i = (float) (j + 1);
1988 dphir += (-phirp_poi_fpar->GetBinContent(j) * TMath::Cos(i * phirp1)
1989 + phirp_poi_fpar->GetBinContent(j + 4) * TMath::Sin(i * phirp1))
1995 for (
int j = 0; j < 4; j++) {
1996 Float_t i = (float) (j + 1);
1997 dphir += (-phirp_poi_fpar->GetBinContent(j) * TMath::Cos(i * phirp2)
1998 + phirp_poi_fpar->GetBinContent(j + 4) * TMath::Sin(i * phirp2))
2003 delrp = (phirp1 - phirp2);
2007 cout <<
"<D-poi> Impact parameter " <<
fMCEventHeader->GetB() <<
", delrp = " << delrp << endl;
2010 delrp = delrp * RADDEG;
2011 if (delrp < -180.) delrp += 360.;
2012 if (delrp > 180.) delrp -= 360.;
2016 phirp = RADDEG * atan2(Qy1 + Qy2, Qx1 + Qx2);
2017 while (phirp < -180.) {
2020 while (phirp > 180.) {
2024 TH1F* phirp_poi_fpar =
fflowFile->Get<TH1F>(
"phirp_poi_fpar");
2026 for (
int j = 0; j < 4; j++) {
2027 Float_t i = (float) (j + 1);
2030 dphir += ((-phirp_poi_fpar->GetBinContent(j) * TMath::Cos(i * phirp / RADDEG)
2031 + phirp_poi_fpar->GetBinContent(j + 4) * TMath::Sin(i * phirp / RADDEG))
2036 phirp += dphir * RADDEG;
2037 while (phirp < -180.) {
2040 while (phirp > 180.) {
2046 while (delrp < -180.) {
2049 while (delrp > 180.) {
2067 Float_t t0m_hit = 0.;
2069 Float_t t0mf_hit = 0.;
2071 Float_t t0mnf_hit = 0.;
2076 Int_t nFTofHits = 0;
2077 Float_t T0FAST[NFHITS];
2078 for (Int_t l = 0; l < NFHITS; l++) {
2082 for (Int_t j = 0; j <
nTofHits; j++) {
2087 if (NULL == TofHit)
continue;
2088 t_hit = TofHit->GetTime();
2093 Float_t dist = TMath::Sqrt(TMath::Power(TofHit->GetX(), 2) + TMath::Power(TofHit->GetY(), 2)
2094 + TMath::Power(TofHit->GetZ(), 2));
2095 Float_t t0_hit = t_hit - dist /
clight;
2096 for (Int_t l = 0; l < NFHITS; l++) {
2097 if (t0_hit < T0FAST[l]) {
2099 if (T0FAST[l] < 100.) {
2100 for (Int_t ll = NFHITS - 1; ll >= l; ll--) {
2101 T0FAST[ll + 1] = T0FAST[ll];
2117 Int_t nfhmax = nFTofHits;
2120 Int_t lmax = NFHITS;
2123 for (Int_t l = 0; l < NFHITS; l++) {
2124 if (T0FAST[l] < 100. && nfh < nfhmax) {
2128 T02 += TMath::Power(T0FAST[l], 2.);
2132 T0MIN = T0MIN / (Float_t) nfh;
2133 Float_t T0RMS = TMath::Sqrt(T02 - T0MIN * T0MIN);
2137 while (T0RMS > 1. && nT0It < 10 && nfh > 2) {
2139 Float_t T0AV = T0MIN;
2140 if (TMath::Abs(T0FAST[lmin] - T0AV) > TMath::Abs(T0FAST[lmax] - T0AV)) {
2150 for (Int_t l = lmin; l < lmax; l++) {
2151 if (T0FAST[l] < 100. && nfh < nfhmax) {
2154 T02 += TMath::Power(T0FAST[l], 2.);
2158 T0MIN = T0MIN / (Float_t) nfh;
2159 T0RMS = TMath::Sqrt(T02 - T0MIN * T0MIN);
2165 for (Int_t j = 0; j <
nTofHits; j++) {
2167 if (NULL == TofHit)
continue;
2168 if (TofHit->GetTime() < 0.2)
continue;
2177 if (NULL == poiMatch) {
2178 LOG(warn) <<
"No MC point found for hit " << j <<
", digi " << iDigInd0;
2181 LOG(debug) <<
"Got PoiMatch for TofHit " << j <<
", digi " << iDigInd0 <<
": " << poiMatch;
2187 LOG(info) <<
"Got invalid PoiMatch for TofHit " << j <<
", digi " << iDigInd0 <<
": " << poiMatch;
2195 LOG(debug) <<
"No Link to MCTofPoint found for hit " << j;
2201 Int_t k = TofPoint->GetTrackID();
2204 px_MC = MCTrack->
GetPx();
2205 py_MC = MCTrack->
GetPy();
2206 pz_MC = MCTrack->
GetPz();
2207 p_MC =
sqrt(px_MC * px_MC + py_MC * py_MC + pz_MC * pz_MC);
2209 cout <<
"<E-hit> Too many MCTracks " << k <<
" from " <<
nMCTracks << endl;
2212 LOG(debug) << Form(
"HadronAnalysis:: hit %d, poi %d, MCt %d, Eff %d, mom %7.2f, pdg %d ", j, lp, k, TrackP[k],
2223 if (TrackP[k] == 0) {
2226 t_hit = TofHit->GetTime();
2227 Float_t delta_tof = TofPoint->GetTime() - t_hit;
2228 Float_t delta_x = TofPoint->GetX() - TofHit->GetX();
2229 Float_t delta_y = TofPoint->GetY() - TofHit->GetY();
2230 Float_t delta_z = TofPoint->GetZ() - TofHit->GetZ();
2232 fa_dxx->Fill(TofPoint->GetX(), delta_x);
2233 fa_dxy->Fill(TofPoint->GetY(), delta_x);
2234 fa_dxz->Fill(TofPoint->GetZ(), delta_x);
2235 fa_dyx->Fill(TofPoint->GetX(), delta_y);
2236 fa_dyy->Fill(TofPoint->GetY(), delta_y);
2237 fa_dyz->Fill(TofPoint->GetZ(), delta_y);
2238 fa_dzx->Fill(TofPoint->GetX(), delta_z);
2239 fa_dzy->Fill(TofPoint->GetY(), delta_z);
2240 fa_dzz->Fill(TofPoint->GetZ(), delta_z);
2242 fa_xy_hit1->Fill(TofHit->GetX(), TofHit->GetY());
2245 fa_dhit_ch->Fill(TofHit->GetCh(), TofHit->GetFlag());
2247 Float_t vel = TofPoint->GetLength() / t_hit;
2248 Float_t bet = vel /
clight;
2249 if (bet > 0.99999) { bet = 0.99999; }
2250 Float_t tofmass = p_MC / bet * TMath::Sqrt(1. - bet * bet) * TMath::Sign(1, pdgCode);
2251 Float_t dist = TMath::Sqrt(TMath::Power(TofHit->GetX(), 2) + TMath::Power(TofHit->GetY(), 2)
2252 + TMath::Power(TofHit->GetZ(), 2));
2257 Float_t t0_hit = t_hit - dist /
clight;
2260 if (t0_hit < T0MIN + 2.4 * MaxT0) {
2262 t0m_hit = ((Float_t)(NT0 - 1) * t0m_hit + t0_hit) / (Float_t) NT0;
2263 if ((TMath::Abs(TofHit->GetX()) < 55. && TMath::Abs(TofHit->GetY()) < 55.)) {
2265 t0mf_hit = ((Float_t)(NT0F - 1) * t0mf_hit + t0_hit) / (Float_t) NT0F;
2269 t0mnf_hit = ((Float_t)(NT0NF - 1) * t0mnf_hit + t0_hit) / (Float_t) NT0NF;
2274 if (TrackP[k] > 1)
continue;
2278 Float_t Phip = RADDEG * atan2(MCTrack->
GetPy(), MCTrack->
GetPx());
2280 if (dphi < -180.) { dphi += 360.; };
2281 if (dphi > 180.) { dphi -= 360.; };
2282 dphi = dphi / RADDEG;
2292 if (use_pions_for_flow && TMath::Abs((MCTrack->
GetRapidity() - yrp_mid) / yrp_mid) >
GetDY()
2311 if (use_pions_for_flow && TMath::Abs((MCTrack->
GetRapidity() - yrp_mid) / yrp_mid) >
GetDY()
2451 if (rp_weight != 0.) {
2452 if (gRandom->Uniform(1) > 0.5) {
2454 Qx1 = Qx1 + rp_weight * MCTrack->
GetPx();
2455 Qy1 = Qy1 + rp_weight * MCTrack->
GetPy();
2459 Qx2 = Qx2 + rp_weight * MCTrack->
GetPx();
2460 Qy2 = Qy2 + rp_weight * MCTrack->
GetPy();
2466 cout <<
"<W> CHA: >=2. hit from track k =" << k <<
" Hit# "
2476 if (Np1 > 0 && Np2 > 0) {
2477 phirp1 = atan2(Qy1, Qx1);
2478 phirp2 = atan2(Qy2, Qx2);
2480 TH1F* phirp_hit_fpar =
fflowFile->Get<TH1F>(
"phirps_hit_fpar");
2482 for (
int jj = 0; jj < 4; jj++) {
2483 Float_t i = (float) (jj + 1);
2484 dphir += (-phirp_hit_fpar->GetBinContent(jj) * TMath::Cos(i * phirp1)
2485 + phirp_hit_fpar->GetBinContent(jj + 4) * TMath::Sin(i * phirp1))
2491 for (
int jj = 0; jj < 4; jj++) {
2492 Float_t i = (float) (jj + 1);
2493 dphir += (-phirp_hit_fpar->GetBinContent(jj) * TMath::Cos(i * phirp2)
2494 + phirp_hit_fpar->GetBinContent(jj + 4) * TMath::Sin(i * phirp2))
2501 delrp = (phirp1 - phirp2);
2504 cout <<
"<D-hit> Impact parameter " <<
fMCEventHeader->GetB() <<
", delrp = " << delrp << endl;
2507 delrp = delrp * RADDEG;
2508 if (delrp < -180.) delrp += 360.;
2509 if (delrp > 180.) delrp -= 360.;
2513 phirp = RADDEG * atan2(Qy1 + Qy2, Qx1 + Qx2);
2514 while (phirp < -180.) {
2517 while (phirp > 180.) {
2521 TH1F* phirp_hit_fpar =
fflowFile->Get<TH1F>(
"phirps_hit_fpar");
2523 for (
int jj = 0; jj < 4; jj++) {
2524 Float_t i = (float) (jj + 1);
2527 dphir += ((-phirp_hit_fpar->GetBinContent(jj) * TMath::Cos(i * phirp / RADDEG)
2528 + phirp_hit_fpar->GetBinContent(jj + 4) * TMath::Sin(i * phirp / RADDEG))
2533 phirp += dphir * RADDEG;
2534 while (phirp < -180.) {
2537 while (phirp > 180.) {
2543 while (delrp < -180.) {
2546 while (delrp > 180.) {
2580 const Double_t DISTMAX = 100.;
2581 Double_t WMAX = 10.;
2584 Int_t NGTofTrack = 0;
2599 cout << Form(
" TofTrack selection: %d. iteration, reassigned hits: %d, "
2600 "GlobTracks %d, TofTracks %d",
2607 if (NRIt == 1) NTHMUL[i] = 0;
2615 cout <<
"<Di> NRIt " << NRIt <<
": Global Track " << i <<
", TofHit " << j <<
" StsTrk " << s << endl;
2619 if (0 == tparf->GetQp()) {
2621 cout <<
"<W> Global Track " << i <<
" without Qp!, take from Sts " << s << endl;
2624 if (NULL == StsTrack) {
2625 cout <<
"<E> Invalid StsTrack pointer at location " << s << endl;
2631 if (0 == tparf->GetQp()) {
2632 cout <<
"<E> Global Track " << i <<
" without Qp! " << endl;
2637 cout <<
"<E> Invalid StsTrack PID " << StsTrack->
GetPidHypo() <<
" at " << s << endl;
2641 FairTrackParam paramExtr;
2648 cout << Form(
"<D> Extrapolate Glob Track %d to prim. vertex %6.2f with "
2656 Float_t momf = 1. / tparf->GetQp();
2657 if (momf < 0.) momf = -momf;
2664 Float_t DISTMIN = WMAX;
2666 Weight_THMUL[i][0] = WMAX;
2671 cout <<
"<Dt> Global Track " << i <<
", TofHit " << j <<
", StsTrk " << s <<
", TofTrack " << tt
2677 if (NULL == TofHit)
continue;
2682 cout <<
"<E> invalid dist for gt " << i <<
", tt " << tt <<
", d:" << dist << endl;
2694 Float_t moml = 1. / tpar->GetQp();
2695 if (moml < 0.) moml = -moml;
2696 Float_t bet = TofHit->GetR() / TofHit->GetTime() /
clight;
2698 if (bet > 0.99999) { bet = 0.99999; }
2699 Float_t tofmass = momf / bet *
sqrt(1. - bet * bet) * TMath::Sign(1., tpar->GetQp());
2701 if (TofTrack->
GetMass() != tofmass) { cout <<
"<E> did not store tofmass properly " << tofmass << endl; }
2703 Float_t mindelmass = 1.E6;
2709 for (Int_t im = 0; im <
NMASS; im++) {
2710 Float_t delmass = TMath::Abs(TMath::Abs(TofTrack->
GetMass()) -
refMass[im]);
2711 if (delmass < mindelmass) {
2712 mindelmass = delmass;
2717 Float_t delp = TMath::Abs((momf - moml) / momf);
2736 cout << Form(
"<D> w for gt %3d, tt %3d, w: %9.5f, d: %7.2f, m: "
2737 "%7.3f, l: %7.2f, dp: %7.3f, p: %7.2f ",
2738 i, tt, w, dist, mindelmass, minlen, delp, momf)
2743 if (nth == MAXNHT) {
2745 cout <<
"<W> Too many TofTrack candidates for track " << i <<
", limit to " << nth << endl;
2751 for (Int_t jth = 0; jth < nth; jth++) {
2753 cout <<
" Compare for position " << jth <<
" w " << w <<
" - " << Weight_THMUL[i][jth] << endl;
2755 if (w < Weight_THMUL[i][jth]) {
2761 cout <<
" Put Track " << tt <<
" with w = " << w <<
" to position " << jthmin <<
" of " << nth
2764 for (Int_t jth = nth; jth > jthmin; jth--) {
2766 cout <<
" Save Track " << IndTHMUL[i][jth - 1] <<
" with w " << Weight_THMUL[i][jth - 1]
2767 <<
" to position " << jth << endl;
2769 IndTHMUL[i][jth] = IndTHMUL[i][jth - 1];
2770 Weight_THMUL[i][jth] = Weight_THMUL[i][jth - 1];
2772 IndTHMUL[i][jthmin] = tt;
2773 Weight_THMUL[i][jthmin] = w;
2777 BestTofTrack = TofTrack;
2780 if (verbose > 5) { cout << Form(
"<DMin> gt %d, hit %d, tt %d, w: %6.2f", i, Bthi, Btt, w) << endl; }
2783 if (verbose > 10) { cout << Form(
"<D> tt-loop: gt %d, tt %d, w: %6.2f", i, tt, w) << endl; }
2786 NTHMUL[i] = nth + 1;
2791 for (Int_t k = 0; k < NTHMUL[i]; k++) {
2793 cout << Form(
"<Ddeb> i %d, k %d, M %d, Ind %d ", i, k, NTHMUL[i], IndTHMUL[i][k]) << endl;
2797 cout <<
"<DSum> GlobTrack " << i <<
", TMul: " << NTHMUL[i] <<
", TofTrack " << IndTHMUL[i][k]
2801 << TofTrack->
GetMass() <<
", mom " << momf <<
", w " << Weight_THMUL[i][k] << endl;
2806 if (NTHMUL[i] > 0) {
2807 Btt = IndTHMUL[i][0];
2808 if (Btt < 0 || Btt >
fTofTracks->GetEntriesFast()) {
2809 cout <<
"<E> invalid TofTrackIndex " << Btt <<
", gt " << i <<
", NRIt " << NRIt << endl;
2817 cout <<
"<DBest> GloBTrack " << i <<
", TMul: " << NTHMUL[i] <<
", TofTrack " << IndTHMUL[i][0]
2818 <<
", TofHit " << BestTofTrack->
GetTofHitIndex() <<
", TMul_hit "
2821 << Weight_THMUL[i][0] << endl;
2830 cout << Form(
"<Ddis> NRIt %d, gt %d, BestTofTrack j=%d, best 0x%p, %d, "
2832 NRIt, i, j, BestTofTrack, Btt, Bthi, Weight_THMUL[i][0])
2836 if (NRIt == 1) GlobTrack->
SetLength(0.);
2840 cout << Form(
"<Ddeb> BestTofTrack j=%d, best 0x%p, %d", j, BestTofTrack, BestTofTrack->
GetTofHitIndex())
2845 cout <<
"<E2> invalid dist for gt " << i <<
", Btt " << Btt <<
", d:" << dist << endl;
2848 if (dist < DISTMAX && Weight_THMUL[i][0] < WMAX) {
2849 Int_t jh = NTofHitTMul[Bthi]++;
2850 Int_t nht = NTofHitTMul[Bthi];
2851 if (nht == MAXNHT) {
2852 if (verbose > -1) { cout <<
"<E> Too many TofTrack candidates for hit " << Bthi <<
", break!" << endl; }
2855 IndTofTrack_TofHit[Bthi][jh] = Btt;
2857 cout <<
"<Ias> GlobTrack " << i <<
" -> TofTrack " << Btt <<
", TofHitIndex " << Bthi <<
", TMul_hit "
2867 cout <<
"<D> GlobTrack " << i <<
": update TofHitIndex from " << j <<
" (Mul " << NTofHitTMul[j] <<
") "
2870 <<
", m " << BestTofTrack->
GetMass() <<
", w " << Weight_THMUL[i][0] <<
", cur: tt "
2871 << IndTofTrack_TofHit[Bthi][0] <<
", gt " << io <<
", w " << Weight_THMUL[io][0] <<
" ? " << endl;
2875 if (Weight_THMUL[i][0] < Weight_THMUL[io][0]) {
2877 cout <<
"<D> New cand. is better, invalidate entry for gt " << io << endl;
2880 NTofHitTMul[Bthi]--;
2881 IndTofTrack_TofHit[Bthi][0] = Btt;
2887 cout << Form(
"<D> Stick to old assignment, Bthi %d, TM %d, THM %d", Bthi, NTofHitTMul[Bthi],
2891 NTofHitTMul[Bthi]--;
2892 if (NTHMUL[i] > 1) {
2894 for (Int_t jth = 0; jth < NTHMUL[i]; jth++) {
2895 IndTHMUL[i][jth] = IndTHMUL[i][jth + 1];
2896 Weight_THMUL[i][jth] = Weight_THMUL[i][jth + 1];
2898 Btt = IndTHMUL[i][0];
2904 cout <<
"<I> no TofTrack candidate for Global Track " << i << endl;
2914 cout <<
"<D> Old choice better, current options: NTHMUL " << NTHMUL[i]
2915 <<
", take btt = " << IndTHMUL[i][0] <<
", bthi " << Bthi << endl;
2921 cout <<
"<E> GlobTrack " << i
2922 <<
": double assignment of hit, check all possibilities "
2936 cout <<
"<D> GlobTrack " << i <<
", dist = " << dist <<
", w = " << Weight_THMUL[i][0]
2937 <<
" -> remove TofTrack" << endl;
2945 cout <<
"<D> GlobTrack " << i <<
", dist = " << dist <<
", w = " << Weight_THMUL[i][0] <<
" -> no TofTrack"
2951 if (verbose > 10) { cout <<
"<Dch> GlobTrack " << i <<
"(" <<
nGlobTracks <<
"), Btt " << Btt << endl; }
2954 cout <<
"<Q> Reassignment iteration for b= " <<
fMCEventHeader->GetB() <<
": " << NReas << endl;
2963 if (j > -1 && Weight_THMUL[i][0] >= WMAX) {
2964 cout <<
"<E> TofHit assigned beyond w-limit, Track " << i <<
" w= " << Weight_THMUL[i][0] << endl;
2968 if (verbose > 10) { cout <<
"<Da> gt " << i <<
", th " << j <<
", s " << s << endl; }
2970 if (0 == tparf->GetQp()) {
2971 if (verbose > 10) cout <<
"<W2> Global Track " << i <<
" without Qp!, take from Sts " << s << endl;
2976 Float_t qpf = tparf->GetQp();
2978 cout <<
"<E2> GlobTrack " << i <<
", STS " << s <<
", TofHit " << j <<
" without momentum " << endl;
2983 Double_t vtxb = 100.;
2986 Int_t NStsMCc[100] = {100 * 0};
2991 FairTrackParam paramExtr;
3000 for (Int_t ih = 0; ih < NStsHits; ih++) {
3003 LOG(debug1) <<
" inspect STS track " << s <<
", hit " << ih <<
", hitindex " << iHind;
3004 if (NULL ==
fStsHits) LOG(fatal) <<
" No STS Hits available ";
3007 if (NULL == hit)
continue;
3008 LOG(debug1) <<
" valid hit " << ih <<
", hitindex " << iHind
3013 std::stringstream ss;
3015 for (Int_t iDigi = 0; iDigi < fclu->
GetNofDigis(); iDigi++) {
3016 ss << fclu->
GetDigi(iDigi) <<
" ";
3021 for (Int_t iL = 0; iL < stsdigiMatch->
GetNofLinks(); iL++) {
3024 if (NULL == poi)
continue;
3025 Int_t MCInd = poi->GetTrackID();
3026 ss <<
" MCInd " << poi->GetTrackID() <<
" ";
3028 for (; iMCt < NStsMCt; iMCt++) {
3029 if (MCInd == StsMCt[iMCt]) {
3034 if (iMCt == NStsMCt) {
3035 LOG(debug) <<
"contribution by new MC track: " << MCInd;
3036 StsMCt[iMCt] = MCInd;
3043 for (Int_t iDigi = 0; iDigi < bclu->
GetNofDigis(); iDigi++) {
3044 ss << bclu->
GetDigi(iDigi) <<
" ";
3049 for (Int_t iL = 0; iL < stsdigiMatch->
GetNofLinks(); iL++) {
3052 Int_t MCInd = poi->GetTrackID();
3053 ss <<
" MCInd " << poi->GetTrackID() <<
" ";
3055 for (; iMCt < NStsMCt; iMCt++) {
3056 if (MCInd == StsMCt[iMCt]) {
3061 if (iMCt == NStsMCt) {
3062 LOG(debug) <<
"contribution by new back MC track: " << MCInd;
3063 StsMCt[iMCt] = MCInd;
3071 LOG(debug1) << ss.str();
3074 std::stringstream ss;
3075 ss <<
"STS summary: NStsMCt =" << NStsMCt;
3076 for (Int_t iT = 0; iT < NStsMCt; iT++) {
3077 ss <<
" iT " << iT <<
" NMCc " << NStsMCc[iT] <<
" MCt " << StsMCt[iT];
3079 LOG(debug) << ss.str();
3083 LOG(debug) <<
"StsTrack " << s <<
" with " << NStsHits <<
" Hits without StsPoints ??? from Global Track " << i
3084 <<
", TofHit " << j;
3088 Int_t iMaxCount = 0;
3089 for (Int_t k = 0; k < NStsMCt; k++) {
3090 if (NStsMCc[k] > iMaxCount) {
3092 iMaxCount = NStsMCc[k];
3103 px_MC = MCTrack->
GetPx();
3104 py_MC = MCTrack->
GetPy();
3105 pz_MC = MCTrack->
GetPz();
3106 p_MC =
sqrt(px_MC * px_MC + py_MC * py_MC + pz_MC * pz_MC);
3172 if (NULL == TofHit)
continue;
3173 if (TofHit->GetTime() < 0.2)
continue;
3178 LOG(fatal) <<
"No Digi Info available for TofHit !?? ";
3188 Int_t iPoiArr[iDigiMul];
3189 Int_t iTrkArr[iDigiMul];
3192 for (Int_t iLink = 0; iLink < digiMatch->
GetNofLinks(); iLink++) {
3196 LOG(error) <<
"MC-Points not available for DigInd " << iDigInd;
3201 LOG(info) <<
"Got PoiMatch for TofHit " << j <<
", t " << TofHit->GetTime() <<
", digi " << iDigInd <<
": "
3203 if (NULL == poiMatch)
continue;
3210 LOG(info) <<
"Got invalid PoiMatch for TofHit " << j <<
", digi " << iDigInd <<
": " << poiMatch;
3215 if (lp != iPoiArr[iPoiMul]) {
3218 iPoiArr[iPoiMul] = lp;
3220 iPoiArr[iPoiMul] = lp;
3224 k = TofPoint->GetTrackID();
3225 if (k != iTrkArr[iTrkMul]) {
3226 iTrkArr[iTrkMul] = k;
3228 iTrkArr[iTrkMul] = k;
3241 LOG(error) <<
"MC-Tracks not available for k " << k;
3245 if (NULL == MCTrack) {
3246 LOG(warn) <<
"No Track for TofHit " << j <<
", Tofpoint " << lp <<
", TrackId " << k;
3250 px_MC = MCTrack->
GetPx();
3251 py_MC = MCTrack->
GetPy();
3252 pz_MC = MCTrack->
GetPz();
3253 p_MC =
sqrt(px_MC * px_MC + py_MC * py_MC + pz_MC * pz_MC);
3258 if (tpar->GetQp() == 0.) {
3259 cout <<
"Invalid momentum for global track " << i <<
" TofHit " << j << endl;
3262 Double_t mom = 1. / tpar->GetQp();
3263 if (mom < 0.) mom = -mom;
3264 Float_t vel = TofHit->GetR() / TofHit->GetTime();
3265 Float_t bet = vel /
clight;
3266 Double_t m2 = mom * mom * (1. / bet / bet - 1.);
3268 if (bet > 0.99999) { bet = 0.99999; }
3269 Float_t tofmass = mom / bet *
sqrt(1. - bet * bet) * TMath::Sign(1., tpar->GetQp());
3281 fa_xy_glo1->Fill(TofHit->GetX(), TofHit->GetY());
3284 fa_m2mom_glo->Fill(mom * TMath::Sign(1., tpar->GetQp()), m2);
3291 if (NULL != TofPoint) {
3313 Float_t Phip = RADDEG * atan2(MCTrack->
GetPy(), MCTrack->
GetPx());
3315 if (dphi < -180.) { dphi += 360.; };
3316 if (dphi > 180.) { dphi -= 360.; };
3317 dphi = dphi / RADDEG;
3332 if (use_pions_for_flow && TMath::Abs((MCTrack->
GetRapidity() - yrp_mid) / yrp_mid) >
GetDY()
3357 if (use_pions_for_flow && TMath::Abs((MCTrack->
GetRapidity() - yrp_mid) / yrp_mid) >
GetDY()
3401 fa_xy_glo_p->Fill(TofHit->GetX(), TofHit->GetY());
3440 fa_xy_glo_d->Fill(TofHit->GetX(), TofHit->GetY());
3461 fa_xy_glo_t->Fill(TofHit->GetX(), TofHit->GetY());
3482 fa_xy_glo_h->Fill(TofHit->GetX(), TofHit->GetY());
3503 fa_xy_glo_a->Fill(TofHit->GetX(), TofHit->GetY());
3541 if (rp_weight != 0.) {
3542 if (gRandom->Uniform(1) > 0.5) {
3544 Qx1 = Qx1 + rp_weight * MCTrack->
GetPx();
3545 Qy1 = Qy1 + rp_weight * MCTrack->
GetPy();
3549 Qx2 = Qx2 + rp_weight * MCTrack->
GetPx();
3550 Qy2 = Qy2 + rp_weight * MCTrack->
GetPy();
3556 if (verbose > 10) { cout <<
"<D> RP analysis " << Np1 <<
", " << Np2 << endl; }
3557 if (Np1 > 0 && Np2 > 0) {
3558 phirp1 = atan2(Qy1, Qx1);
3559 phirp2 = atan2(Qy2, Qx2);
3561 TH1F* phirp_glo_fpar =
fflowFile->Get<TH1F>(
"phirps_glo_fpar");
3563 for (
int j = 0; j < 4; j++) {
3564 Float_t i = (float) (j + 1);
3565 dphir += (-phirp_glo_fpar->GetBinContent(j) * TMath::Cos(i * phirp1)
3566 + phirp_glo_fpar->GetBinContent(j + 4) * TMath::Sin(i * phirp1))
3572 for (
int j = 0; j < 4; j++) {
3573 Float_t i = (float) (j + 1);
3574 dphir += (-phirp_glo_fpar->GetBinContent(j) * TMath::Cos(i * phirp2)
3575 + phirp_glo_fpar->GetBinContent(j + 4) * TMath::Sin(i * phirp2))
3580 delrp = (phirp1 - phirp2);
3586 delrp = delrp * RADDEG;
3587 if (delrp < -180.) delrp += 360.;
3588 if (delrp > 180.) delrp -= 360.;
3591 phirp = RADDEG * atan2(Qy1 + Qy2, Qx1 + Qx2);
3592 while (phirp < -180.) {
3595 while (phirp > 180.) {
3599 TH1F* phirp_glo_fpar =
fflowFile->Get<TH1F>(
"phirp_glo_fpar");
3601 for (
int j = 0; j < 4; j++) {
3602 Float_t i = (float) (j + 1);
3605 dphir += ((-phirp_glo_fpar->GetBinContent(j) * TMath::Cos(i * phirp / RADDEG)
3606 + phirp_glo_fpar->GetBinContent(j + 4) * TMath::Sin(i * phirp / RADDEG))
3611 phirp += dphir * RADDEG;
3612 while (phirp < -180.) {
3615 while (phirp > 180.) {
3621 while (delrp < -180.) {
3624 while (delrp > 180.) {
3638 LOG(info) <<
"-I- CbmHadronAnalysis::Exec : "
3639 <<
"event " <<
fEvents <<
" processed.";
3697#include "TLorentzVector.h"
3698#include "TVector3.h"
3700 static TH1F* fhTofChi;
3701 static TH1F* fhDperp;
3702 static TH2F* fhdEdxMul;
3703 static TH2F* fhdEdxMulsec;
3704 static TH1F* fhDTRDprim;
3705 static TH1F* fhDTRDsec;
3706 static TH1F* fhDperp2;
3707 static TH1F* fhDperpS;
3708 static TH1F* fhD0prim;
3709 static TH1F* fhOpAng;
3711 static TH1F* fhMinv;
3712 static TH1F* fhPathLen;
3713 static TH1F* fhMMom;
3714 static TH1F* fhMIXOpAng;
3715 static TH1F* fhMIXDCA;
3716 static TH1F* fhMIXMinv;
3717 static TH1F* fhMIXPathLen;
3718 static TH1F* fhMIXMMom;
3719 static TH1F* fhMCLamMom;
3720 static TH1F* fhMCPathLen;
3722 static TH2F* fa_ptm_rap_gen_lam;
3723 static TH2F* fa_ptm_rap_rec_lam;
3724 static TH2F* fa_ptm_rap_mix_lam;
3727 static TH2F* fhTofStsXX;
3728 static TH2F* fhTofStsYY;
3729 static TH2F* fhTofStsdXdY;
3730 static TH2F* fhTofStsdXX;
3731 static TH2F* fhTofStsdXY;
3732 static TH2F* fhTofStsdYX;
3733 static TH2F* fhTofStsdYY;
3734 static TH2F* fhTofVelNhit;
3736 static TH2F* fhStsStsX0Y0;
3737 static TH2F* fhStsStsX0X1;
3738 static TH2F* fhStsStsY0X1;
3739 static TH2F* fhStsStsX0Y1;
3740 static TH2F* fhStsStsY0Y1;
3741 static TH2F* fhStsStsX0X2;
3742 static TH2F* fhStsStsY0X2;
3743 static TH2F* fhStsStsX0Y2;
3744 static TH2F* fhStsStsY0Y2;
3746 static TH3F* fhTofSts1dXXY;
3747 static TH3F* fhTofSts1dYXY;
3748 static TH3F* fhTofSts2dXXY;
3749 static TH3F* fhTofSts2dYXY;
3751 static std::vector<std::list<std::vector<TLorentzVector>>> fvP;
3752 static std::vector<std::list<std::vector<TLorentzVector>>> fvX;
3753 static std::vector<std::list<std::vector<TVector3>>> fvX0;
3754 static std::vector<std::list<std::vector<TVector3>>> fvDX;
3769 const Int_t fiNMixClasses = 10;
3770 const Double_t beamRotY = -25.;
3771 const Double_t MLAM = 1.1156;
3772 const Double_t DMLAM = 0.015;
3773 const Int_t NSECMASS = 2;
3774 const Int_t iMode = 0;
3775 Float_t secMass[NSECMASS] = {0.139, 0.938};
3780 secMass[1] = 2.808381;
3783 const Double_t dTofSigX = 0.5;
3784 const Double_t dTofSigY = 0.8;
3785 const Double_t dTofSigT = 0.08;
3786 const Double_t dChiTofLim = 3.;
3793 LOG(debug) <<
"Secondaries from " <<
nTofHits <<
" TofHits, " <<
nStsHits <<
" StsHits and " << nTrdHits
3794 <<
" TrdHits in event " <<
iCandEv;
3798 TDirectory* oldDir = gDirectory;
3799 fHist->ReOpen(
"Update");
3801 Double_t MinvMin = secMass[0] + secMass[1];
3802 cout <<
"Add secondary histos to " <<
fHist->GetName() << endl;
3803 fhTofChi =
new TH1F(Form(
"hTofChi"), Form(
"TofHit Merger; #chi "), 100, 0., dChiTofLim * 2.);
3804 fhDperp =
new TH1F(Form(
"hDperp"), Form(
"transverse matching distance; d [cm]"), 100, 0., 5.);
3806 new TH2F(Form(
"hdEdxMul"), Form(
"average energy loss vs TrdHitMul; TrdHitMul; dE []"), 4, 1, 5, 100, 0., 5.E-5);
3807 fhdEdxMulsec =
new TH2F(Form(
"hdEdxMulsec"), Form(
"av. energy loss vs TrdHitMul for secondaries; TrdHitMul; dE []"),
3808 4, 1, 5, 100, 0., 5.E-5);
3810 new TH1F(Form(
"hDTRDprim"), Form(
"TRD transverse matching distance (prim); d [cm]"), 100, 0., 2. *
fdDistTRD);
3812 new TH1F(Form(
"hDTRDsec"), Form(
"TRD transverse matching distance (sec); d [cm]"), 100, 0., 2. *
fdDistTRD);
3814 fhDperp2 =
new TH1F(Form(
"hDperp2"), Form(
"transverse matching distance (prim); d [cm]"), 100, 0., 1.);
3815 fhDperpS =
new TH1F(Form(
"hDperpS"), Form(
"transverse matching distance (sec); d [cm]"), 100, 0., 1.);
3816 fhD0prim =
new TH1F(Form(
"hD0prim"), Form(
"transverse distance to primary vertex; d [cm]"), 100, 0., 2.);
3817 fhOpAng =
new TH1F(Form(
"hOpAng"), Form(
"opening angle; #alpha [rad]"), 100, 0., 0.5);
3818 fhDCA =
new TH1F(Form(
"hDCA"), Form(
"distance of closest approach; d [cm]"), 100, 0., 2.);
3819 fhMinv =
new TH1F(Form(
"hMinv"), Form(
"invariant mass; M_{inv} [GeV]"), 100, MinvMin, MinvMin + 0.2);
3820 fhPathLen =
new TH1F(Form(
"hPathLen"), Form(
"path length; L [cm]"), 100, 0., 30.);
3821 fhMMom =
new TH1F(Form(
"hMMom"), Form(
"momentum of mother ; p [GeV]"), 100, 0., 5.);
3822 fhMIXOpAng =
new TH1F(Form(
"hMIXOpAng"), Form(
"opening angle; #alpha [rad]"), 100, 0., 0.5);
3823 fhMIXDCA =
new TH1F(Form(
"hMIXDCA"), Form(
"distance of closest approach; d [cm]"), 100, 0., 2.);
3824 fhMIXMinv =
new TH1F(Form(
"hMIXMinv"), Form(
"invariant mass; M_{inv} [GeV]"), 100, MinvMin, MinvMin + 0.2);
3825 fhMIXPathLen =
new TH1F(Form(
"hMIXPathLen"), Form(
"path length; L [cm]"), 100, 0., 30.);
3826 fhMIXMMom =
new TH1F(Form(
"hMIXMMom"), Form(
"momentum of mother ; p [GeV]"), 100, 0., 5.);
3827 fhMCPathLen =
new TH1F(Form(
"hMCPathLen"), Form(
"MC hyperon path length; L [cm]"), 100, 0., 30.);
3828 fhMCLamMom =
new TH1F(Form(
"hMCLamMom"), Form(
"MC hyperon momentum; p [GeV]"), 100, 0., 5.);
3830 fhTofStsXX =
new TH2F(Form(
"hTofStsXX"), Form(
"Position correlation XX; x_{Sts} [cm]; x_{TofExt} [cm]"), 200, -10.,
3831 10., 200, -15., 15.);
3832 fhTofStsYY =
new TH2F(Form(
"hTofStsYY"), Form(
"Position correlation YY; y_{Sts} [cm]; y_{TofExt} [cm]"), 200, -10.,
3833 10., 200, -15., 15.);
3835 new TH2F(Form(
"hTofStsdXdY"), Form(
"Position correlation dXX; #Deltax_{Tof-Sts} [cm]; #Deltay_{Tof-Sts} [cm]"),
3836 100, -2., 2., 100, -2., 2.);
3837 fhTofStsdXX =
new TH2F(Form(
"hTofStsdXX"), Form(
"Position correlation dXX; x_{Tof} [cm]; #Deltax_{Tof-Sts} [cm]"),
3838 200, -15., 15., 100, -2., 2.);
3839 fhTofStsdXY =
new TH2F(Form(
"hTofStsdXY"), Form(
"Position correlation dXY; y_{Tof} [cm]; #Deltax_{Tof-Sts} [cm]"),
3840 200, -15., 15., 100, -2., 2.);
3841 fhTofStsdYX =
new TH2F(Form(
"hTofStsdYX"), Form(
"Position correlation dYX; x_{Tof} [cm]; #Deltay_{Tof-Sts} [cm]"),
3842 200, -15., 15., 100, -2., 2.);
3843 fhTofStsdYY =
new TH2F(Form(
"hTofStsdYY"), Form(
"Position correlation dYY; y_{Tof} [cm]; #Deltay_{Tof-Sts} [cm]"),
3844 200, -15., 15., 100, -2., 2.);
3845 fhTofVelNhit =
new TH2F(Form(
"hTofVelNhit"), Form(
"Velocity vs Nhit; Nhit []; v [cm/ns]"), 10, 0, 10, 100, 0, 50.);
3846 fhNT0 =
new TH1F(Form(
"hNT0"), Form(
"Number of T0 signal per event; N []; cts []"), 10, 0, 10);
3847 fhStsStsX0Y0 =
new TH2F(Form(
"hStsStsX0Y0"), Form(
"vertex location; x0_{StsSts} [cm]; y0_{StsSts} [cm]"), 100, -2.,
3850 fhStsStsX0X1 =
new TH2F(Form(
"hStsStsX0X1"), Form(
"vertex position dependence; x1_{Sts} [cm]; x0_{StsSts} [cm]"),
3851 100, -10., 10., 100, -2., 2.);
3852 fhStsStsY0X1 =
new TH2F(Form(
"hStsStsY0X1"), Form(
"vertex position dependence; x1_{Sts} [cm]; y0_{StsSts} [cm]"),
3853 100, -10., 10., 100, -2., 2.);
3854 fhStsStsX0Y1 =
new TH2F(Form(
"hStsStsX0Y1"), Form(
"vertex position dependence; y1_{Sts} [cm]; x0_{StsSts} [cm]"),
3855 100, -10., 10., 100, -2., 2.);
3856 fhStsStsY0Y1 =
new TH2F(Form(
"hStsStsY0Y1"), Form(
"vertex position dependence; y1_{Sts} [cm]; y0_{StsSts} [cm]"),
3857 100, -10., 10., 100, -2., 2.);
3859 fhStsStsX0X2 =
new TH2F(Form(
"hStsStsX0X2"), Form(
"vertex position dependence; x2_{Sts} [cm]; x0_{StsSts} [cm]"),
3860 100, -10., 10., 100, -2., 2.);
3861 fhStsStsY0X2 =
new TH2F(Form(
"hStsStsY0X2"), Form(
"vertex position dependence; x2_{Sts} [cm]; y0_{StsSts} [cm]"),
3862 100, -10., 10., 100, -2., 2.);
3863 fhStsStsX0Y2 =
new TH2F(Form(
"hStsStsX0Y2"), Form(
"vertex position dependence; y2_{Sts} [cm]; x0_{StsSts} [cm]"),
3864 100, -10., 10., 100, -2., 2.);
3865 fhStsStsY0Y2 =
new TH2F(Form(
"hStsStsY0Y2"), Form(
"vertex position dependence; y2_{Sts} [cm]; y0_{StsSts} [cm]"),
3866 100, -10., 10., 100, -2., 2.);
3869 new TH3F(Form(
"hTofSts1dXXY"), Form(
"Position deviation dXXY; x_{Sts1} [cm]; y_{Sts1}; #Deltax_{Tof-Sts} [cm]"),
3870 100, -10., 10., 100, 10, 10, 50, -2., 2.);
3872 new TH3F(Form(
"hTofSts1dYXY"), Form(
"Position deviation dYXY; x_{Sts1} [cm]; y_{Sts1}; #Deltay_{Tof-Sts} [cm]"),
3873 100, -10., 10., 100, 10, 10, 50, -2., 2.);
3875 new TH3F(Form(
"hTofSts2dXXY"), Form(
"Position deviation dXXY; x_{Sts1} [cm]; y_{Sts1}; #Deltax_{Tof-Sts} [cm]"),
3876 100, -10., 10., 100, 10, 10, 50, -2., 2.);
3878 new TH3F(Form(
"hTofSts2dYXY"), Form(
"Position deviation dYXY; x_{Sts1} [cm]; y_{Sts1}; #Deltay_{Tof-Sts} [cm]"),
3879 100, -10., 10., 100, 10, 10, 50, -2., 2.);
3884 Float_t ptmmax = 2.5;
3888 fa_ptm_rap_gen_lam =
3889 new TH2F(
"ptm_rap_gen_lam",
"MCTrack-gen lam; y; p_{T}/m", ptm_nbx, ymin, ymax, ptm_nby, 0., ptmmax);
3890 fa_ptm_rap_rec_lam =
new TH2F(
"ptm_rap_rec_lam",
"rec lam; y; p_{T}/m", ptm_nbx, ymin, ymax, ptm_nby, 0., ptmmax);
3891 fa_ptm_rap_mix_lam =
new TH2F(
"ptm_rap_mix_lam",
"mix lam; y; p_{T}/m", ptm_nbx, ymin, ymax, ptm_nby, 0., ptmmax);
3893 gDirectory = oldDir;
3902 if (MCTrack->
GetMotherId() == -1 && pdgCode == 3122) {
3903 fhMCLamMom->Fill(MCTrack->
GetP());
3905 PLAM.RotateY(beamRotY * TMath::Pi() / 180.);
3906 fa_ptm_rap_gen_lam->Fill(PLAM.Rapidity(), PLAM.Pt() / MCTrack->
GetMass());
3908 if (MCTrack->
GetMotherId() > -1 && pdgCode == -211) {
3913 fhMCPathLen->Fill(MCV.Mag());
3914 LOG(debug) <<
"MC vertex at Pathlen = " << MCV.Mag() <<
", pi-mom " << MCTrack->
GetP();
3919 fvP.resize(fiNMixClasses);
3920 fvX.resize(fiNMixClasses);
3921 fvX0.resize(fiNMixClasses);
3922 fvDX.resize(fiNMixClasses);
3925 if (iMixClass >= fiNMixClasses) iMixClass = fiNMixClasses - 1;
3927 std::vector<TLorentzVector> P;
3929 std::vector<TLorentzVector> X;
3931 std::vector<TVector3> X0;
3933 std::vector<TVector3> DX;
3935 const Int_t NTrdStations = 4;
3936 std::vector<std::vector<Int_t>> iTRD;
3946 Double_t dTrdDistMin[
nTofHits][NTrdStations];
3948 Int_t proton_cand = 0;
3949 Int_t pion_cand = 0;
3951 for (Int_t j = 0; j <
nStsHits; j++) {
3953 dTofDistMin[j] = 100.;
3954 dTofDist2Min[j] = 100.;
3960 for (Int_t i = 0; i <
nTofHits; i++) {
3962 if (NULL == pTofHit)
continue;
3963 if (pTofHit->
GetZ() == 0.) {
3971 LOG(debug) <<
"Found T0 " << dT0 <<
" from " << dNT0 <<
" signals ";
3977 for (Int_t i = 0; i <
nTofHits; i++) {
3979 if (NULL == pTofHit)
continue;
3980 if (pTofHit->
GetZ() == 0.)
continue;
3986 for (Int_t i = 0; i <
nTofHits; i++) {
3988 if (NULL == pTofHit)
continue;
3989 if (pTofHit->
GetZ() == 0.)
continue;
3990 for (Int_t i2 = 0; i2 <
nTofHits; i2++) {
3993 if (NULL == pTofHit2)
continue;
3994 if (pTofHit2->
GetZ() == 0.)
continue;
3996 if (pTofHit2->
GetZ() < pTofHit->
GetZ()) {
3999 pTofHit2 = pTofHittmp;
4001 Double_t dPosZExp = pTofHit->
GetZ() / pTofHit2->
GetZ();
4002 Double_t dPosXExp = pTofHit2->
GetX() * dPosZExp;
4003 Double_t dPosYExp = pTofHit2->
GetY() * dPosZExp;
4004 Double_t dTimeExp = pTofHit2->
GetTime() * dPosZExp;
4005 Double_t dChi2 = TMath::Power((dPosXExp - pTofHit->
GetX()) / dTofSigX, 2)
4006 + TMath::Power((dPosYExp - pTofHit->
GetY()) / dTofSigY, 2)
4007 + TMath::Power((dTimeExp - pTofHit->
GetTime()) / dTofSigT, 2);
4008 Double_t dChi = TMath::Sqrt(dChi2) / 3.;
4009 fhTofChi->Fill(dChi);
4010 if (dChi < dChiTofLim) {
4017 LOG(debug) <<
"Tof Hits " << i <<
" and " << i2 <<
" merged ";
4018 LOG(debug) <<
"Tof " << i <<
", xyz " << pTofHit->
GetX() <<
", " << pTofHit->
GetY() <<
", "
4028 for (Int_t i = 0; i <
nTofHits; i++) {
4030 if (NULL == pTofHit) {
4031 P[i].SetPxPyPzE(0., 0., 0., 0.);
4034 if (pTofHit->
GetZ() == 0.)
continue;
4035 dStsDistMin[i] = 1.E3;
4036 dSts2DistMin[i] = 1.E3;
4037 Double_t dTofX = 0, dTofY = 0, dStsX = 0, dStsY = 0;
4038 for (Int_t l = 0; l < NTrdStations; l++)
4039 dTrdDistMin[i][l] = 1.E3;
4042 for (Int_t j = 0; j <
nStsHits; j++) {
4045 Double_t sPosZ = pStsHit->
GetZ();
4046 Double_t sPosXext = pTofHit->
GetX() * sPosZ / pTofHit->
GetZ();
4047 Double_t sPosYext = pTofHit->
GetY() * sPosZ / pTofHit->
GetZ();
4048 Double_t dDist2 = TMath::Power(pStsHit->
GetX() - sPosXext, 2) + TMath::Power(pStsHit->
GetY() - sPosYext, 2);
4049 Double_t dDist = TMath::Sqrt(dDist2);
4051 LOG(debug) <<
"Sts " << j <<
", xyz " << pStsHit->
GetX() <<
", " << pStsHit->
GetY() <<
", " << sPosZ;
4052 LOG(debug) <<
"Tof " << i <<
", xyz " << pTofHit->
GetX() <<
", " << pTofHit->
GetY() <<
", " << pTofHit->
GetZ();
4053 LOG(debug) <<
"Tof " << i <<
", Sts " << j << Form(
" -> dist %6.3f, Min %6.3f ", dDist, dStsDistMin[i]);
4056 if (iTofMin[j] > -1) {
4057 LOG(debug) << Form(
"Sts hit %d already assigned to tof hit %d with "
4058 "dist= %6.3f, prev %6.3f",
4059 j, iTofMin[j], dDist, dTofDistMin[j]);
4060 if (dDist > dTofDistMin[j])
continue;
4062 dStsDistMin[i] = dDist;
4064 dTofDistMin[j] = dDist;
4065 LOG(debug) <<
"Prim Track cand started for Tof " << i <<
", Sts " << j
4066 << Form(
": dist %6.3f, Min %6.3f at z = %4.1f", dDist, dStsDistMin[i], sPosZ);
4069 dStsX = pStsHit->
GetX();
4070 dStsY = pStsHit->
GetY();
4074 if (dStsDistMin[i] < 100.) {
4075 fhDperp->Fill(dStsDistMin[i]);
4076 fhTofStsXX->Fill(dStsX, dTofX);
4077 fhTofStsYY->Fill(dStsY, dTofY);
4078 fhTofStsdXdY->Fill(dTofX - dStsX, dTofY - dStsY);
4079 fhTofStsdXX->Fill(dTofX, dTofX - dStsX);
4080 fhTofStsdXY->Fill(dTofY, dTofX - dStsX);
4081 fhTofStsdYX->Fill(dTofX, dTofY - dStsY);
4082 fhTofStsdYY->Fill(dTofY, dTofY - dStsY);
4087 for (Int_t j = 0; j <
nStsHits; j++) {
4088 if (iTofMin[j] < 0)
continue;
4089 Int_t i = iTofMin[j];
4092 for (Int_t k = 0; k <
nStsHits; k++) {
4093 if (j == k)
continue;
4095 Double_t sPos2Z = pSts2Hit->
GetZ();
4096 Double_t sPos2Xext = pStsHit->
GetX() * sPos2Z / pStsHit->
GetZ();
4097 Double_t sPos2Yext = pStsHit->
GetY() * sPos2Z / pStsHit->
GetZ();
4098 Double_t dDist2 = TMath::Power(pSts2Hit->
GetX() - sPos2Xext, 2) + TMath::Power(pSts2Hit->
GetY() - sPos2Yext, 2);
4099 Double_t dDist = TMath::Sqrt(dDist2);
4100 LOG(debug) <<
"Tof " << i <<
", Sts " << j
4101 << Form(
" Sts2 %d -> dist %6.3f, Min %6.3f at z = %4.1f", k, dDist, dSts2DistMin[i], sPos2Z);
4104 if (iTofMin[k] > -1) {
4105 LOG(debug) << Form(
"Sts2hit %d already assigned to tof hit %d with "
4106 "dist= %6.3f, prev %6.3f",
4107 k, iTofMin[k], dDist, dTofDistMin[k]);
4108 if (dDist > dTofDist2Min[k])
continue;
4110 dSts2DistMin[i] = dDist;
4112 if (pStsHit->
GetZ() < pSts2Hit->
GetZ()) {
4120 dTofDistMin[k] = dDist;
4121 LOG(debug) <<
"Prim Track cand found for Tof " << i <<
", Sts " << j <<
", Sts2 " << k
4122 << Form(
": dist %6.3f, Min %6.3f at z = %4.1f", dDist, dSts2DistMin[i], sPos2Z);
4126 fhDperp2->Fill(dSts2DistMin[i]);
4132 fhStsStsX0Y0->Fill(dX0, dY0);
4133 fhStsStsX0X1->Fill(pSts0->
GetX(), dX0);
4134 fhStsStsX0Y1->Fill(pSts0->
GetY(), dX0);
4135 fhStsStsY0X1->Fill(pSts0->
GetX(), dY0);
4136 fhStsStsY0Y1->Fill(pSts0->
GetY(), dY0);
4137 fhStsStsX0X2->Fill(pSts1->
GetX(), dX0);
4138 fhStsStsX0Y2->Fill(pSts1->
GetY(), dX0);
4139 fhStsStsY0X2->Fill(pSts1->
GetX(), dY0);
4140 fhStsStsY0Y2->Fill(pSts1->
GetY(), dY0);
4143 fhTofSts1dXXY->Fill(pSts0->
GetX(), pSts0->
GetY(),
4145 fhTofSts1dYXY->Fill(pSts0->
GetX(), pSts0->
GetY(),
4147 fhTofSts2dXXY->Fill(pSts1->
GetX(), pSts1->
GetY(),
4149 fhTofSts2dYXY->Fill(pSts1->
GetX(), pSts1->
GetY(),
4156 for (Int_t i = 0; i <
nTofHits; i++) {
4157 Int_t j = iStsMin[i][1];
4158 if (j < 0)
continue;
4160 if (NULL == pTofHit)
continue;
4161 if (pTofHit->
GetZ() == 0)
continue;
4163 Double_t dDx = pTofHit->
GetX() - pStsHit->
GetX();
4164 Double_t dDy = pTofHit->
GetY() - pStsHit->
GetY();
4165 Double_t dDz = pTofHit->
GetZ() - pStsHit->
GetZ();
4166 LOG(debug) <<
"Check for TRD hits between STS " << j <<
" and TOF " << i;
4168 for (Int_t l = 0; l < nTrdHits; l++) {
4171 Double_t dXexp = pStsHit->
GetX() + dDx * (pTrdHit->
GetZ() - pStsHit->
GetZ()) / dDz;
4172 Double_t dYexp = pStsHit->
GetY() + dDy * (pTrdHit->
GetZ() - pStsHit->
GetZ()) / dDz;
4174 TMath::Sqrt(TMath::Power(dXexp - pTrdHit->
GetX(), 2) + TMath::Power(dYexp - pTrdHit->
GetY(), 2));
4176 LOG(debug) <<
"Inspect TRD hit " << l <<
" in "
4177 << Form(
"Module 0x%08x, layer %d", pTrdHit->
GetAddress(),
4179 <<
" at z= " << pTrdHit->
GetZ() <<
" dD = " << dDtrans <<
" < " <<
fdDistTRD;
4180 fhDTRDprim->Fill(dDtrans);
4181 if (dDtrans <
fdDistTRD && dDtrans < dTrdDistMin[i][iTrdLayer]) {
4182 Int_t iMul = iTRD[i].size();
4183 if (dTrdDistMin[i][iTrdLayer] < 1.E3) {
4186 for (; ll < iMul; ll++)
4191 dTrdDistMin[i][iTrdLayer] = dDtrans;
4192 iTRD[i].resize(iMul + 1);
4195 LOG(debug) <<
"assign TrdHit " << l <<
" to TofHit " << i <<
" in layer " << iTrdLayer
4196 <<
" with d = " << dDtrans <<
", TrdMul" << iMul <<
", dEdx = " << pTrdHit->
GetELoss();
4201 for (Int_t i = 0; i <
nTofHits; i++) {
4202 Int_t iMul = iTRD[i].size();
4204 Double_t ddEdx = 0.;
4205 for (Int_t l = 0; l < iMul; l++) {
4209 ddEdx /= (Double_t) iMul;
4210 fhdEdxMul->Fill((Double_t) iMul, ddEdx);
4216 for (Int_t i = 0; i <
nTofHits; i++) {
4218 if (NULL == pTofHit)
continue;
4219 if (pTofHit->
GetZ() == 0.)
continue;
4220 if (iStsMin[i][0] > -1 && iStsMin[i][1] > -1) {
4224 Double_t dDx = pStsHit->
GetX() - pSts2Hit->
GetX();
4225 Double_t dDy = pStsHit->
GetY() - pSts2Hit->
GetY();
4226 Double_t dDz = pStsHit->
GetZ() - pSts2Hit->
GetZ();
4228 Double_t dX0 = pSts2Hit->
GetX() - dDx / dDz * pSts2Hit->
GetZ();
4229 Double_t dY0 = pSts2Hit->
GetY() - dDy / dDz * pSts2Hit->
GetZ();
4230 Double_t dD0 = TMath::Sqrt(dX0 * dX0 + dY0 * dY0);
4231 fhD0prim->Fill(dD0);
4235 Double_t dDd = TMath::Sqrt(dDx * dDx + dDy * dDy + dDz * dDz);
4236 Double_t vel = pTofHit->
GetR() / pTofHit->
GetTime();
4237 Double_t bet = vel /
clight;
4238 if (bet > 0.9999)
continue;
4239 Double_t m = secMass[1];
4240 Double_t pmag = m * bet / TMath::Sqrt(1. - bet * bet);
4241 Double_t pz = pmag * dDz / dDd;
4242 Double_t px = pmag * dDx / dDd;
4243 Double_t py = pmag * dDy / dDd;
4244 Double_t E = TMath::Sqrt(pmag * pmag + m * m);
4245 P[i].SetPxPyPzE(px, py, pz, E);
4247 LOG(debug) <<
"Init proton LV at ind " << Form(
"%d %d %d", i, iStsMin[i][0], iStsMin[i][1])
4248 <<
" with beta = " << bet <<
", minv = " << P[i].M() <<
", tof " << X[i].T() <<
", TRDHmul "
4250 X0[i].SetXYZ(pSts2Hit->
GetX(), pSts2Hit->
GetY(), pSts2Hit->
GetZ());
4251 DX[i].SetXYZ(dDx, dDy, dDz);
4258 for (Int_t i = 0; i <
nTofHits; i++) {
4260 if (NULL == pTofHit)
continue;
4261 if (pTofHit->
GetZ() == 0.)
continue;
4262 LOG(debug) <<
"Tof " << i << Form(
" sec cand Min %6.3f > %6.3f ?", dStsDistMin[i],
fdDistPrimLim);
4264 Double_t dDistMin = 100.;
4267 for (Int_t j = 0; j <
nStsHits; j++) {
4268 LOG(debug) <<
"Tof " << i <<
", Sts " << j
4269 << Form(
" ? sec cand %6.3f Min %6.3f ", dTofDistMin[j],
fdDistPrimLim);
4273 Double_t dDx = pTofHit->
GetX() - pStsHit->
GetX();
4274 Double_t dDy = pTofHit->
GetY() - pStsHit->
GetY();
4275 Double_t dDz = pTofHit->
GetZ() - pStsHit->
GetZ();
4277 for (Int_t k = 0; k <
nStsHits; k++) {
4278 if (j == k)
continue;
4280 Double_t sPos2Z = pSts2Hit->
GetZ();
4281 Double_t sPos2Xext = pStsHit->
GetX() + dDx / dDz * (sPos2Z - pStsHit->
GetZ());
4282 Double_t sPos2Yext = pStsHit->
GetY() + dDy / dDz * (sPos2Z - pStsHit->
GetZ());
4284 TMath::Power(pSts2Hit->
GetX() - sPos2Xext, 2) + TMath::Power(pSts2Hit->
GetY() - sPos2Yext, 2);
4285 Double_t dDist = TMath::Sqrt(dDist2);
4286 LOG(debug) <<
"Sec Tof " << i <<
", Sts " << j
4287 << Form(
" Sts2 %d -> dist %6.3f < %6.3f ? at z = %4.1f", k, dDist,
4297 fhDperpS->Fill(dDistMin);
4298 LOG(debug) <<
"Sec Dist for TofHit " << i <<
": " << dDistMin <<
", j " << jbest <<
", k " << kbest;
4300 if (dDistMin < 100.) {
4303 if (pSts2Hit->
GetZ() > pStsHit->
GetZ()) {
4309 Double_t dDx = pTofHit->
GetX() - pStsHit->
GetX();
4310 Double_t dDy = pTofHit->
GetY() - pStsHit->
GetY();
4311 Double_t dDz = pTofHit->
GetZ() - pStsHit->
GetZ();
4312 for (Int_t l = 0; l < nTrdHits; l++) {
4315 Double_t dXexp = pStsHit->
GetX() + dDx * (pTrdHit->
GetZ() - pStsHit->
GetZ()) / dDz;
4316 Double_t dYexp = pStsHit->
GetY() + dDy * (pTrdHit->
GetZ() - pStsHit->
GetZ()) / dDz;
4318 TMath::Sqrt(TMath::Power(dXexp - pTrdHit->
GetX(), 2) + TMath::Power(dYexp - pTrdHit->
GetY(), 2));
4320 fhDTRDsec->Fill(dDtrans);
4321 LOG(debug) <<
"Inspect sec. TRD hit " << l <<
" in "
4322 << Form(
"Module 0x%08x, layer %d", pTrdHit->
GetAddress(),
4324 <<
" at z= " << pTrdHit->
GetZ() <<
" dD = " << dDtrans <<
" < " <<
fdDistTRD;
4326 && dDtrans < dTrdDistMin[i][iTrdLayer]) {
4327 Int_t iMul = iTRD[i].size();
4328 if (dTrdDistMin[i][iTrdLayer] < 1.E3) {
4331 for (; ll < iMul; ll++)
4337 dTrdDistMin[i][iTrdLayer] = dDtrans;
4338 iTRD[i].resize(iMul + 1);
4341 LOG(debug) <<
"assign TrdHit " << l <<
" to TofHit " << i <<
" in layer " << iTrdLayer
4342 <<
" with d = " << dDtrans <<
", TrdMul" << iMul <<
", dEdx = " << pTrdHit->
GetELoss();
4346 Int_t iMul = iTRD[i].size();
4348 Double_t ddEdx = 0.;
4349 for (Int_t l = 0; l < iMul; l++) {
4353 ddEdx /= (Double_t) iMul;
4354 fhdEdxMulsec->Fill((Double_t) iMul, ddEdx);
4358 Double_t dDx = pStsHit->
GetX() - pSts2Hit->
GetX();
4359 Double_t dDy = pStsHit->
GetY() - pSts2Hit->
GetY();
4360 Double_t dDz = pStsHit->
GetZ() - pSts2Hit->
GetZ();
4361 Double_t dDd = TMath::Sqrt(dDx * dDx + dDy * dDy + dDz * dDz);
4362 Double_t vel = pTofHit->
GetR() / pTofHit->
GetTime();
4363 Double_t bet = vel /
clight;
4364 if (bet > 0.9999)
continue;
4365 Double_t m = secMass[0];
4366 Double_t pmag = m * bet / TMath::Sqrt(1. - bet * bet);
4367 Double_t pz = pmag * dDz / dDd;
4368 Double_t px = pmag * dDx / dDd;
4369 Double_t py = pmag * dDy / dDd;
4370 Double_t E = TMath::Sqrt(pmag * pmag + m * m);
4371 P[i].SetPxPyPzE(px, py, pz, E);
4373 LOG(debug) <<
"Init pion LV at ind " << i <<
" with beta = " << bet <<
", minv = " << P[i].M() <<
", tof "
4374 << X[i].T() <<
", TRDHmul " << iTRD[i].size();
4375 X0[i].SetXYZ(pSts2Hit->
GetX(), pSts2Hit->
GetY(), pSts2Hit->
GetZ());
4376 DX[i].SetXYZ(dDx, dDy, dDz);
4382 LOG(debug) <<
" Ev " <<
iCandEv <<
" has " << proton_cand <<
" protons and " << pion_cand <<
" pion candidates";
4383 if (proton_cand > 0 && pion_cand > 0) {
4384 LOG(debug) <<
"add event " <<
iCandEv <<
" to mixing class " << iMixClass <<
" of size " << fvP[iMixClass].size();
4386 fvP[iMixClass].push_front(P);
4387 fvX[iMixClass].push_front(X);
4388 fvX0[iMixClass].push_front(X0);
4389 fvDX[iMixClass].push_front(DX);
4392 fvP[iMixClass].pop_back();
4393 fvX[iMixClass].pop_back();
4394 fvX0[iMixClass].pop_back();
4395 fvDX[iMixClass].pop_back();
4402 for (Int_t i = 0; i <
nTofHits; i++) {
4403 if (TMath::Abs(P[i].M() - secMass[1]) < 0.01) {
4404 std::list<std::vector<TLorentzVector>>::iterator itX = fvX[iMixClass].begin();
4405 std::list<std::vector<TVector3>>::iterator itX0 = fvX0[iMixClass].begin();
4406 std::list<std::vector<TVector3>>::iterator itDX = fvDX[iMixClass].begin();
4411 LOG(debug) <<
"LV P has size " << P.size() <<
", fvP size " << fvP[iMixClass].size() <<
" in mix class "
4413 for (std::list<std::vector<TLorentzVector>>::iterator itP = fvP[iMixClass].begin(); itP != fvP[iMixClass].end();
4419 std::vector<TLorentzVector> PE = *itP;
4420 std::vector<TLorentzVector> XE = *itX;
4421 std::vector<TVector3> X0E = *itX0;
4422 std::vector<TVector3> DXE = *itDX;
4423 LOG(debug) <<
"iMixEv " << iMixEv <<
": PE has size " << PE.size() <<
", X0E: " << X0E.size();
4426 for (UInt_t j = 0; j < PE.size(); j++) {
4427 LOG(debug) <<
"P comp " << j <<
": " << P[j].M() <<
" " << PE[j].M();
4429 LOG(debug) <<
"P not properly restored from list ";
4433 for (UInt_t j = 0; j < PE.size(); j++) {
4434 if (TMath::Abs(PE[j].M() - secMass[0]) < 0.01) {
4436 Double_t dOpAngle = DX[i].Angle(DXE[j]);
4437 if (iMixEv == 1) fhOpAng->Fill(dOpAngle);
4439 fhMIXOpAng->Fill(dOpAngle);
4442 TVector3 N = DX[i].Cross(DXE[j]);
4443 if (N.Mag() == 0.)
continue;
4445 Double_t dDCA = TMath::Abs((X0[i] - X0E[j]) * N);
4446 if (iMixEv == 1) fhDCA->Fill(dDCA);
4448 fhMIXDCA->Fill(dDCA);
4449 LOG(debug) <<
"DCA for iMixEv " << iMixEv <<
" at ind i " << i <<
", j " << j <<
": " << dDCA;
4450 if (dDCA == 0.)
continue;
4453 TVector3 D = dDCA * N;
4454 TVector3 Ni = DX[i].Cross(N);
4455 Double_t cj = -(X0E[j] - X0[i] - D) * Ni / (DXE[j] * Ni);
4456 TVector3 V = X0E[j] + cj * DXE[j] - 0.5 * D;
4457 Double_t dVLen = V.Mag();
4460 TLorentzVector PM = P[i] + PE[j];
4461 TVector3 PV = TVector3(PM.Px(), PM.Py(), PM.Pz());
4462 Double_t TofM = dVLen / PM.Beta() /
clight;
4464 Double_t TofMLast = 100.;
4466 std::vector<TLorentzVector> Ptmp = P;
4467 std::vector<TLorentzVector> PEtmp = PE;
4468 std::vector<TLorentzVector> Xtmp = X;
4469 std::vector<TLorentzVector> XEtmp = XE;
4471 while (TMath::Abs(TofM - TofMLast) > 0.001 && Niter++ < 3) {
4472 LOG(debug) <<
"MinvI at ind i " << i <<
", j " << j <<
": " << PM.M() <<
", vertex: " << V[0] <<
" "
4473 << V[1] <<
" " << V[2] <<
", Len " << dVLen <<
", mom = " << PV.Mag() <<
", TofM " << TofM;
4476 for (Int_t ii = 0; ii < 2; ii++) {
4483 TofX = Xtmp[k].Vect();
4489 TofX = XEtmp[k].Vect();
4493 TVector3 vDTofV = TofX - V;
4494 Double_t vel = vDTofV.Mag() / (tof - TofM);
4495 Double_t bet = vel /
clight;
4496 if (bet > 0.9999) bet = 0.9999;
4497 Double_t pmag = m * bet / TMath::Sqrt(1. - bet * bet);
4498 TVector3 vPsec = vDTofV;
4500 Double_t E = TMath::Sqrt(pmag * pmag + m * m);
4502 Ptmp[k].SetVect(vPsec);
4506 PEtmp[k].SetVect(vPsec);
4510 PM = Ptmp[i] + PEtmp[j];
4511 PV = TVector3(PM.Px(), PM.Py(), PM.Pz());
4513 TofM = dVLen / PM.Beta() /
clight;
4515 Double_t minv = PM.M();
4518 fhPathLen->Fill(dVLen);
4519 fhMMom->Fill(PV.Mag());
4521 LOG(debug) <<
"MinvII in event " <<
fEvents <<
" at ind i " << i <<
", j " << j <<
": " << minv
4522 <<
", vertex: " << V[0] <<
" " << V[1] <<
" " << V[2] <<
", Len " << dVLen
4523 <<
", mom = " << PV.Mag() <<
", tof " << TofM;
4525 if (TMath::Abs(MLAM - minv) < DMLAM) {
4526 PM.RotateY(beamRotY * TMath::Pi() / 180.);
4527 fa_ptm_rap_rec_lam->Fill(PM.Rapidity(), PM.Pt() / MLAM);
4531 fhMIXMinv->Fill(minv);
4532 fhMIXPathLen->Fill(dVLen);
4533 fhMIXMMom->Fill(PV.Mag());
4534 if (TMath::Abs(MLAM - minv) < DMLAM) {
4535 PM.RotateY(beamRotY * TMath::Pi() / 180.);
4536 fa_ptm_rap_mix_lam->Fill(PM.Rapidity(), PM.Pt() / MLAM);