30 #include "FairGeoNode.h"
31 #include "FairGeoTransform.h"
32 #include "FairRootManager.h"
33 #include "FairRunAna.h"
34 #include "FairRuntimeDb.h"
38 #include "TClonesArray.h"
47 #include "TMCProcess.h"
64 : FairTask(
"RichGeoTestQa")
70 , fRichRingMatches(NULL)
85 , fDrawEventDisplay(true) {
101 cout <<
"CbmRichGeoTest::Init" << endl;
102 FairRootManager* ioman = FairRootManager::Instance();
104 Fatal(
"CbmRichGeoTest::Init",
"RootManager not instantised!");
109 if (mcManager ==
nullptr)
110 LOG(fatal) <<
"CbmRichGeoTest::Init() NULL MCDataManager.";
113 if (NULL ==
fMcTracks) { LOG(fatal) <<
"CbmRichGeoTest::Init No MCTrack!"; }
117 LOG(fatal) <<
"CbmRichGeoTest::Init No RichPoint!";
122 LOG(fatal) <<
"CbmRichGeoTest::Init No RefPlanePoint!";
125 fRichHits = (TClonesArray*) ioman->GetObject(
"RichHit");
126 if (NULL ==
fRichHits) { LOG(fatal) <<
"CbmRichGeoTest::Init No RichHit!"; }
131 LOG(fatal) <<
"CbmRichGeoTest::Init No RichDigi!";
135 LOG(fatal) <<
"CbmRichGeoTest::Init No RichDigiMatch!";
138 fRichRings = (TClonesArray*) ioman->GetObject(
"RichRing");
139 if (NULL ==
fRichRings) { LOG(fatal) <<
"CbmRichGeoTest::Init No RichRing!"; }
143 LOG(fatal) <<
"CbmRichGeoTest::Init No RichRingMatch!";
148 LOG(fatal) <<
"CbmRichGeoTest::Init No MCEventList!";
161 cout <<
"CbmRichGeoTest, event No. " <<
fEventNum << endl;
186 "fhHitsXY;X [cm];Y [cm];Counter",
194 "fhPointsXY;X [cm];Y [cm];Counter",
202 "fhPointsXYNoRotation;X [cm];Y [cm];Counter",
209 fHM->
Create1<TH1D>(
"fhHitsZ",
"fhHitsZ;Z [cm];Yield", 1000, 150, 250);
210 fHM->
Create1<TH1D>(
"fhPointsZ",
"fhPointsZ;Z [cm];Yield", 100, 190, 250);
212 for (Int_t
i = 0;
i < 2;
i++) {
214 if (
i == 0) ss <<
"_hits";
215 if (
i == 1) ss <<
"_points";
219 "fhNofHits" + t +
";Nof hits in ring;Yield",
225 "fhNofHits" + t +
";Nof points in ring;Yield",
231 "fhBoverAVsMom" + t +
";p [GeV/c];B/A;Yield",
239 "fhXcYcEllipse" + t +
";X [cm];Y [cm];Yield",
247 "fhBaxisVsMom" + t +
";p [GeV/c];B axis [cm];Yield",
255 "fhAaxisVsMom" + t +
";p [GeV/c];A axis [cm];Yield",
263 "fhChi2EllipseVsMom" + t +
";p [GeV/c];#Chi^{2};Yield",
272 "fhXcYcCircle" + t +
";x [cm];y [cm];Yield",
280 "fhRadiusVsMom" + t +
";p [GeV/c];Radius [cm];Yield",
288 "fhChi2CircleVsMom" + t +
";p [GeV/c];#Chi^{2};Yield",
296 "fhDRVsMom" + t +
";p [GeV/c];dR [cm];Yield",
305 "fhBaxisUpHalf" + t +
";B axis [cm];Yield",
310 "fhBaxisDownHalf" + t +
";B axis [cm];Yield",
317 "fhNofPhotonsPerHit;Number of photons per hit;Yield",
325 "fhDiffAaxis;Nof hits in ring;A_{point}-A_{hit} [cm];Yield",
334 "fhDiffBaxis;Nof hits in ring;B_{point}-B_{hit} [cm];Yield",
343 "fhDiffXcEllipse;Nof hits in ring;Xc_{point}-Xc_{hit} [cm];Yield",
352 "fhDiffYcEllipse;Nof hits in ring;Yc_{point}-Yc_{hit} [cm];Yield",
361 "fhDiffXcCircle;Nof hits in ring;Xc_{point}-Xc_{hit} [cm];Yield",
370 "fhDiffYcCircle;Nof hits in ring;Yc_{point}-Yc_{hit} [cm];Yield",
379 "fhDiffRadius;Nof hits in ring;Radius_{point}-Radius_{hit} [cm];Yield",
389 "fhRadiusVsNofHits;Nof hits in ring;Radius [cm];Yield",
397 "fhAaxisVsNofHits;Nof hits in ring;A axis [cm];Yield",
405 "fhBaxisVsNofHits;Nof hits in ring;B axis [cm];Yield",
413 "fhDRVsNofHits;Nof hits in ring;dR [cm];Yield",
423 "fhDiffXhit",
"fhDiffXhit;X_{point}-X_{hit} [cm];Yield", 200, -.5, .5);
425 "fhDiffYhit",
"fhDiffYhit;Y_{point}-Y_{hit} [cm];Yield", 200, -.5, .5);
429 "fhNofHitsAll",
"fhNofHitsAll;Nof hits in ring;Yield", 50, 0., 50.);
431 "fhNofHitsCircleFit;Nof hits in ring;Yield",
436 "fhNofHitsEllipseFit;Nof hits in ring;Yield",
441 "fhNofHitsCircleFitEff;Nof hits in ring;Efficiency [%]",
446 "fhNofHitsEllipseFitEff;Nof hits in ring;Efficiency [%]",
452 fHM->
Create1<TH1D>(
"fhMcMomEl",
"fhMcMomEl;p [GeV/c];Yield", 24, 0., 12.);
454 "fhMcPtyEl;Rapidity;P_{t} [GeV/c];Yield",
461 fHM->
Create1<TH1D>(
"fhAccMomEl",
"fhAccMomEl;p [GeV/c];Yield", 24, 0., 12.);
463 "fhAccPtyEl;Rapidity;P_{t} [GeV/c];Yield",
471 fHM->
Create1<TH1D>(
"fhMcMomPi",
"fhMcMomPi;p [GeV/c];Yield", 24, 0., 12.);
473 "fhMcPtyPi;Rapidity;P_{t} [GeV/c];Yield",
480 fHM->
Create1<TH1D>(
"fhAccMomPi",
"fhAccMomPi;p [GeV/c];Yield", 24, 0., 12.);
482 "fhAccPtyPi;Rapidity;P_{t} [GeV/c];Yield",
492 "fhNofHitsXYZ;X [cm];Y [cm];Nof hits in ring",
503 "fhNofPointsXYZ;X [cm];Y [cm];Nof points in ring",
514 "fhBoverAXYZ;X [cm];Y [cm];B/A",
525 "fhBaxisXYZ;X [cm];Y [cm];B axis [cm]",
536 "fhAaxisXYZ;X [cm];Y [cm];A axis [cm]",
547 "fhRadiusXYZ;X [cm];Y [cm];Radius [cm]",
558 "fhdRXYZ;X [cm];Y [cm];dR [cm]",
577 "fhNofHitsVsX;X [cm];Nof hits in ring",
585 "fhNofHitsVsY;Abs(Y) [cm];Nof hits in ring",
594 "fhNofPointsVsX;X [cm];Nof points in ring",
602 "fhNofPointsVsY;Abs(Y) [cm];Nof points in ring",
611 "fhBoverAVsX;X [cm];B/A",
619 "fhBoverAVsY;Abs(Y) [cm];B/A",
628 "fhBaxisVsX;X [cm];B axis [cm]",
636 "fhBaxisVsY;Abs(Y) [cm];B axis [cm]",
645 "fhAaxisVsX;X [cm];A axis [cm]",
653 "fhAaxisVsY;Abs(Y) [cm];A axis [cm]",
662 "fhRadiusVsX;X [cm];Radius [cm]",
670 "fhRadiusVsY;Abs(Y) [cm];Radius [cm]",
679 "fhdRVsX",
"fhdRVsX;X [cm];dR [cm]", nBinsX1, xMin1, xMax1, 100, -1., 1.);
681 "fhdRVsY;Abs(Y) [cm];dR [cm]",
690 "fhPhotonEPlaneZ+;Photon energy [eV];Counter",
695 "fhPhotonEPlaneZ-;Photon energy [eV];Counter",
700 "fhPhotonEPmtPoint;Photon energy [eV];Counter",
705 "fhPhotonEPmtHit;Photon energy [eV];Counter",
711 "fhLambdaPlaneZ+;Photon wavelength [nm];Counter",
716 "fhLambdaPlaneZ-;Photon wavelength [nm];Counter",
721 "fhLambdaPmtPoint;Photon wavelength [nm];Counter",
726 "fhLambdaPmtHit;Photon wavelength [nm];Counter",
734 for (Int_t iE = 0; iE < nofEvents; iE++) {
739 for (Int_t iT = 0; iT < nofMcTracks; iT++) {
742 if (!mcTrack)
continue;
744 Int_t pdg = TMath::Abs(mcTrack->
GetPdgCode());
745 Bool_t isMcPrimaryElectron =
746 (pdg == 11 && motherId == -1)
749 if (isMcPrimaryElectron) {
750 fHM->
H1(
"fhMcMomEl")->Fill(mcTrack->
GetP());
754 if (pdg == 211 && motherId == -1) {
755 fHM->
H1(
"fhMcMomPi")->Fill(mcTrack->
GetP());
761 for (Int_t iP = 0; iP < nofPoints; iP++) {
764 if (point == NULL)
continue;
765 TVector3 inPos(point->GetX(), point->GetY(), point->GetZ());
768 fHM->
H2(
"fhPointsXYNoRotation")->Fill(point->GetX(), point->GetY());
769 fHM->
H2(
"fhPointsXY")->Fill(outPos.X(), outPos.Y());
770 fHM->
H1(
"fhPointsZ")->Fill(point->GetZ());
773 point->Momentum(mom);
774 Double_t momMag = mom.Mag();
775 fHM->
H1(
"fhPhotonEPmtPoint")->Fill(1.e9 * momMag);
777 fHM->
H1(
"fhLambdaPmtPoint")->Fill(lambda);
781 for (Int_t iP = 0; iP < nofPlanePoints; iP++) {
784 if (point == NULL)
continue;
786 point->Momentum(mom);
787 Double_t momMag = mom.Mag();
789 if (point->GetPz() > 0) {
790 fHM->
H1(
"fhPhotonEPlaneZ+")->Fill(1.e9 * momMag);
791 fHM->
H1(
"fhLambdaPlaneZ+")->Fill(lambda);
793 fHM->
H1(
"fhPhotonEPlaneZ-")->Fill(1.e9 * momMag);
794 fHM->
H1(
"fhLambdaPlaneZ-")->Fill(lambda);
802 Int_t nofRings =
fRichRings->GetEntriesFast();
803 for (Int_t iR = 0; iR < nofRings; iR++) {
805 if (NULL == ring)
continue;
808 if (NULL == ringMatch)
continue;
813 if (NULL == mcTrack)
continue;
816 Int_t pdg = TMath::Abs(mcTrack->
GetPdgCode());
817 Double_t momentum = mcTrack->
GetP();
818 Double_t pt = mcTrack->
GetPt();
820 Bool_t isMcPrimaryElectron =
821 (pdg == 11 && motherId == -1)
825 if (isMcPrimaryElectron) {
826 fHM->
H1(
"fhAccMomEl")->Fill(momentum);
827 fHM->
H2(
"fhAccPtyEl")->Fill(rapidity, pt);
829 if (pdg == 211 && motherId == -1) {
830 fHM->
H1(
"fhAccMomPi")->Fill(momentum);
831 fHM->
H2(
"fhAccPtyPi")->Fill(rapidity, pt);
835 if (!isMcPrimaryElectron)
continue;
839 for (
int iPoint = 0; iPoint < nofRichPoints; iPoint++) {
842 if (NULL == richPoint)
continue;
843 Int_t trackId = richPoint->GetTrackID();
844 if (trackId < 0)
continue;
847 if (NULL == mcTrackRich)
continue;
849 if (motherIdRich == mcTrackId) {
851 richPoint->Position(posPoint);
891 fHM->
H3(
"fhNofHitsXYZ")->Fill(xc, yc, nh);
892 fHM->
H3(
"fhNofPointsXYZ")->Fill(xc, yc, np);
893 fHM->
H3(
"fhBoverAXYZ")->Fill(xc, yc, b / a);
894 fHM->
H3(
"fhBaxisXYZ")->Fill(xc, yc, b);
895 fHM->
H3(
"fhAaxisXYZ")->Fill(xc, yc, a);
896 fHM->
H3(
"fhRadiusXYZ")->Fill(xc, yc, r);
898 fHM->
H2(
"fhNofHitsVsX")->Fill(xc, nh);
899 fHM->
H2(
"fhNofPointsVsX")->Fill(xc, np);
900 fHM->
H2(
"fhBoverAVsX")->Fill(xc, b / a);
901 fHM->
H2(
"fhBaxisVsX")->Fill(xc, b);
902 fHM->
H2(
"fhAaxisVsX")->Fill(xc, a);
903 fHM->
H2(
"fhRadiusVsX")->Fill(xc, r);
905 fHM->
H2(
"fhNofHitsVsY")->Fill(abs(yc), nh);
906 fHM->
H2(
"fhNofPointsVsY")->Fill(abs(yc), np);
907 fHM->
H2(
"fhBoverAVsY")->Fill(abs(yc), b / a);
908 fHM->
H2(
"fhBaxisVsY")->Fill(abs(yc), b);
909 fHM->
H2(
"fhAaxisVsY")->Fill(abs(yc), a);
910 fHM->
H2(
"fhRadiusVsY")->Fill(abs(yc), r);
912 for (
int iH = 0; iH < ringHit.
GetNofHits(); iH++) {
915 double dr = r -
sqrt((xc - xh) * (xc - xh) + (yc - yh) * (yc - yh));
916 fHM->
H3(
"fhdRXYZ")->Fill(xc, yc, dr);
917 fHM->
H2(
"fhdRVsX")->Fill(xc, dr);
918 fHM->
H2(
"fhdRVsY")->Fill(abs(yc), dr);
934 if (histIndex == 0) {
936 }
else if (histIndex == 1) {
941 fHM->
H1(
"fhBoverAVsMom" + t)->Fill(momentum, axisB / axisA);
942 fHM->
H2(
"fhXcYcEllipse" + t)->Fill(xcEllipse, ycEllipse);
944 fHM->
H1(
"fhNofHits" + t)->Fill(nofHitsRing);
946 if (ycEllipse > 149 || ycEllipse < -149) {
947 fHM->
H1(
"fhBaxisUpHalf" + t)->Fill(axisB);
949 fHM->
H1(
"fhBaxisDownHalf" + t)->Fill(axisB);
952 fHM->
H2(
"fhBaxisVsMom" + t)->Fill(momentum, axisB);
953 fHM->
H2(
"fhAaxisVsMom" + t)->Fill(momentum, axisA);
954 fHM->
H2(
"fhChi2EllipseVsMom" + t)
956 if (histIndex == 0) {
957 fHM->
H2(
"fhAaxisVsNofHits")->Fill(nofHitsRing, axisA);
958 fHM->
H2(
"fhBaxisVsNofHits")->Fill(nofHitsRing, axisB);
971 if (histIndex == 0) {
973 }
else if (histIndex == 1) {
976 fHM->
H1(
"fhXcYcCircle" + t)->Fill(xcCircle, ycCircle);
977 fHM->
H1(
"fhRadiusVsMom" + t)->Fill(momentum, radius);
978 fHM->
H1(
"fhChi2CircleVsMom" + t)
981 for (
int iH = 0; iH < nofHitsRing; iH++) {
985 -
sqrt((xcCircle - xh) * (xcCircle - xh)
986 + (ycCircle - yh) * (ycCircle - yh));
987 fHM->
H1(
"fhDRVsMom" + t)->Fill(momentum, dr);
989 if (histIndex == 0) {
990 fHM->
H2(
"fhDRVsNofHits")->Fill(nofHitsRing, dr);
994 if (histIndex == 0) {
995 fHM->
H2(
"fhRadiusVsNofHits")->Fill(nofHitsRing, radius);
1001 fHM->
H2(
"fhDiffAaxis")
1003 fHM->
H2(
"fhDiffBaxis")
1005 fHM->
H2(
"fhDiffXcEllipse")
1007 fHM->
H2(
"fhDiffYcEllipse")
1013 fHM->
H2(
"fhDiffXcCircle")
1015 fHM->
H2(
"fhDiffYcCircle")
1017 fHM->
H2(
"fhDiffRadius")
1023 Int_t nofHits =
fRichHits->GetEntriesFast();
1024 for (Int_t iH = 0; iH < nofHits; iH++) {
1026 if (hit == NULL)
continue;
1029 if (digiIndex < 0)
continue;
1031 if (NULL == digi)
continue;
1034 if (NULL == digiMatch)
continue;
1036 vector<CbmLink> links = digiMatch->
GetLinks();
1037 for (UInt_t
i = 0;
i < links.size();
i++) {
1038 Int_t pointId = links[
i].GetIndex();
1039 Int_t eventId = links[
i].GetEntry();
1040 if (pointId < 0)
continue;
1044 if (NULL == pMCpt)
continue;
1046 TVector3 inPos(pMCpt->GetX(), pMCpt->GetY(), pMCpt->GetZ());
1049 fHM->
H1(
"fhDiffXhit")->Fill(hit->
GetX() - outPos.X());
1050 fHM->
H1(
"fhDiffYhit")->Fill(hit->
GetY() - outPos.Y());
1053 pMCpt->Momentum(mom);
1054 Double_t momMag = mom.Mag();
1055 fHM->
H1(
"fhPhotonEPmtHit")->Fill(1.e9 * momMag);
1057 fHM->
H1(
"fhLambdaPmtHit")->Fill(lambda);
1071 TCanvas* c =
fHM->
CreateCanvas(ss.str().c_str(), ss.str().c_str(), 500, 500);
1072 c->SetGrid(
true,
true);
1073 TH2D* pad =
new TH2D(ss.str().c_str(),
1074 (ss.str() +
";X [cm];Y [cm]").c_str(),
1081 pad->SetStats(
false);
1086 double xmin = 99999., xmax = -99999., ymin = 99999., ymax = -99999.;
1090 if (xmin > hitX) xmin = hitX;
1091 if (xmax < hitX) xmax = hitX;
1092 if (ymin > hitY) ymin = hitY;
1093 if (ymax < hitY) ymax = hitY;
1095 double xCur = (xmin + xmax) / 2.;
1096 double yCur = (ymin + ymax) / 2.;
1099 TEllipse* circle =
new TEllipse(ringHit->
GetCenterX() - xCur,
1102 circle->SetFillStyle(0);
1103 circle->SetLineWidth(3);
1105 TEllipse* center =
new TEllipse(
1111 TEllipse* hitDr =
new TEllipse(
1113 hitDr->SetFillColor(kRed);
1119 TEllipse* pointDr =
new TEllipse(
1121 pointDr->SetFillColor(kBlue);
1127 ss2 <<
"(r, n)=(" << setprecision(3) << ringHit->
GetRadius() <<
", "
1129 TLatex* latex =
new TLatex(-8., 8., ss2.str().c_str());
1134 Int_t nofMc = (Int_t)
fHM->
H1(
"fhMcMomEl")->GetEntries();
1136 (TH1D*)
fHM->
H1(
"fhNofHits_hits")->Clone(
"fhAccVsMinNofHitsHist");
1137 hist->GetXaxis()->SetTitle(
"Required min nof hits in ring");
1138 hist->GetYaxis()->SetTitle(
"Detector acceptance [%]");
1140 for (
int i = hist->GetNbinsX();
i > 1;
i--) {
1141 sum +=
fHM->
H1(
"fhNofHits_hits")->GetBinContent(
i);
1142 hist->SetBinContent(
i, 100. * sum / nofMc);
1155 (TH1D*) hist->ProjectionY((
string(hist->GetName()) +
"_py").c_str())
1158 py->Scale(1. / py->Integral());
1159 py->GetYaxis()->SetTitle(
"Yield");
1173 fHM->
CreateCanvas(
"richgeo_points_xy",
"richgeo_points_xy", 1200, 1200);
1179 "richgeo_points_xy_no_rotation",
1187 fHM->
H1(
"fhHitsZ")->Scale(1. /
fHM->
H1(
"fhHitsZ")->Integral());
1193 fHM->
H1(
"fhPointsZ")->Scale(1. /
fHM->
H1(
"fhPointsZ")->Integral());
1197 for (
int i = 0;
i < 2;
i++) {
1201 }
else if (
i == 1) {
1206 "richgeo" + t +
"_ellipse_boa_vs_mom");
1207 fHM->
H2(
"fhBoverAVsMom" + t)->GetYaxis()->SetRangeUser(0.8, 1.0);
1211 (
"richgeo" + t +
"_ellipse_xc_yc").c_str(),
1218 "richgeo" + t +
"_chi2_ellipse_vs_mom");
1220 "richgeo" + t +
"_a_vs_mom");
1222 "richgeo" + t +
"_b_vs_mom");
1227 (
"richgeo" + t +
"_b_up_down_halves").c_str(),
1233 (TH1D*)
fHM->
H1(
"fhBaxisUpHalf" + t)->Clone(),
true,
true, 3., 6.);
1236 (TH1D*)
fHM->
H1(
"fhBaxisDownHalf" + t)->Clone(),
true,
true, 3., 6.);
1241 (
"richgeo" + t +
"_circle").c_str(),
1247 cout <<
"Number of hits/points = " <<
fHM->
H1(
"fhNofHits" + t)->GetMean()
1255 "richgeo" + t +
"_chi2_circle_vs_mom");
1257 "richgeo" + t +
"_r_vs_mom");
1259 "richgeo" + t +
"_dr_vs_mom");
1260 fHM->
H2(
"fhDRVsMom" + t)->GetYaxis()->SetRangeUser(-1.05, 1.05);
1265 "richgeo_nof_photons_per_hit",
"richgeo_nof_photons_per_hit", 800, 800);
1266 fHM->
H1(
"fhNofPhotonsPerHit")
1267 ->Scale(1. /
fHM->
H1(
"fhNofPhotonsPerHit")->Integral());
1273 "richgeo_diff_ellipse",
"richgeo_diff_ellipse", 1200, 600);
1295 fHM->
CreateCanvas(
"richgeo_diff_circle",
"richgeo_diff_circle", 900, 600);
1313 fHM->
CreateCanvas(
"richgeo_hits_residual",
"richgeo_hits", 1200, 600);
1316 fHM->
H1(
"fhDiffXhit")->Scale(1. /
fHM->
H1(
"fhDiffXhit")->Integral());
1319 fHM->
H1(
"fhDiffYhit")->Scale(1. /
fHM->
H1(
"fhDiffYhit")->Integral());
1329 (TH1D*)
fHM->
H1(
"fhNofHitsCircleFit")->Clone(),
1330 (TH1D*)
fHM->
H1(
"fhNofHitsEllipseFit")->Clone()},
1331 {
"All",
"Circle fit",
"Ellipse fit"},
1340 TH1D* fhNofHitsCircleFitEff =
1342 TH1D* fhNofHitsEllipseFitEff =
1345 DrawH1(fhNofHitsCircleFitEff);
1346 TLatex* circleFitEffTxt =
new TLatex(
1351 cout <<
"Circle fit efficiency:" << circleFitEffTxt <<
"%" << endl;
1352 circleFitEffTxt->Draw();
1354 DrawH1(fhNofHitsEllipseFitEff);
1355 TLatex* ellipseFitEffTxt =
new TLatex(
1360 cout <<
"Ellipse fit efficiency:" << ellipseFitEffTxt <<
"%" << endl;
1361 ellipseFitEffTxt->Draw();
1386 (TH1D*)
fHM->
H1(
"fhMcMomEl")->Clone(),
1389 "Geometrical acceptance [%]");
1391 (TH2D*)
fHM->
H1(
"fhMcPtyEl")->Clone(),
1394 "Geometrical acceptance [%]");
1397 "richgeo_acc_eff_el_mom",
"richgeo_acc_eff_el_mom", 800, 800);
1399 (TH1D*)
fHM->
H1(
"fhMcMomEl")->Clone());
1400 cout <<
"Geometrical acceptance electrons:" << effEl <<
"%" << endl;
1402 {
"e^{#pm} (" + effEl +
"%)"},
1414 "richgeo_acc_eff_el_pty",
"richgeo_acc_eff_el_pty", 800, 800);
1440 (TH1D*)
fHM->
H1(
"fhMcMomPi")->Clone(),
1443 "Geometrical acceptance [%]");
1445 (TH2D*)
fHM->
H1(
"fhMcPtyPi")->Clone(),
1448 "Geometrical acceptance [%]");
1451 "richgeo_acc_eff_pi_mom",
"richgeo_acc_eff_pi_mom", 800, 800);
1453 (TH1D*)
fHM->
H1(
"fhMcMomPi")->Clone());
1454 cout <<
"Geometrical acceptance pions:" << effPi <<
"%" << endl;
1456 {
"#pi^{#pm} (" + effPi +
"%)"},
1468 "richgeo_acc_eff_pi_pty",
"richgeo_acc_eff_pi_pty", 800, 800);
1470 pyzPiEff->GetZaxis()->SetRangeUser(0, 100);
1475 "richgeo_acc_eff_el_zoom",
"richgeo_acc_eff_el_zoom", 1000, 500);
1478 TH1D* fhMcMomElClone = (TH1D*)
fHM->
H1(
"fhMcMomEl")->Clone();
1479 TH1D* fhAccMomElClone = (TH1D*)
fHM->
H1(
"fhAccMomEl")->Clone();
1480 fhMcMomElClone->GetXaxis()->SetRangeUser(0., 3.);
1481 fhAccMomElClone->GetXaxis()->SetRangeUser(0., 3.);
1482 fhMcMomElClone->SetMinimum(0.);
1483 DrawH1({fhMcMomElClone, fhAccMomElClone},
1492 gPad->SetLogy(
false);
1494 TH1D* px_eff_clone = (TH1D*) pxEff->Clone();
1495 px_eff_clone->GetXaxis()->SetRangeUser(0., 3.);
1496 px_eff_clone->SetMinimum(0.);
1503 "richgeo_numbers_vs_xy_hits",
"richgeo_numbers_vs_xy_hits", 1800, 600);
1515 "richgeo_numbers_vs_xy_points",
1529 "richgeo_numbers_vs_xy_boa",
"richgeo_numbers_vs_xy_boa", 1800, 600);
1535 fHM->
H2(
"fhBoverAVsX")->GetYaxis()->SetRangeUser(0.75, 1.0);
1538 fHM->
H2(
"fhBoverAVsY")->GetYaxis()->SetRangeUser(0.75, 1.0);
1543 "richgeo_numbers_vs_xy_b",
"richgeo_numbers_vs_xy_b", 1800, 600);
1555 "richgeo_numbers_vs_xy_a",
"richgeo_numbers_vs_xy_a", 1800, 600);
1567 "richgeo_numbers_vs_xy_r",
"richgeo_numbers_vs_xy_r", 1800, 600);
1579 "richgeo_numbers_vs_xy_dr",
"richgeo_numbers_vs_xy_dr", 1800, 600);
1591 "richgeo_acc_vs_min_nof_hits",
"richgeo_acc_vs_min_nof_hits", 600, 600);
1593 h->GetXaxis()->SetRangeUser(0., 40.0);
1601 fHM->
H2(
"fhDRVsNofHits")->GetYaxis()->SetRangeUser(-1.05, 1.05);
1609 fHM->
H2(
"fhRadiusVsNofHits")
1611 (
string(
fHM->
H2(
"fhRadiusVsNofHits")->GetName()) +
"_py").c_str()),
1618 fHM->
H2(
"fhAaxisVsNofHits")
1620 (
string(
fHM->
H2(
"fhAaxisVsNofHits")->GetName()) +
"_py").c_str()),
1627 fHM->
H2(
"fhBaxisVsNofHits")
1629 (
string(
fHM->
H2(
"fhBaxisVsNofHits")->GetName()) +
"_py").c_str()),
1639 "richgeo_photon_energy",
"richgeo_photon_energy", 1500, 750);
1643 fHM->
H1(
"fhPhotonEPlaneZ-"),
1644 fHM->
H1(
"fhPhotonEPmtPoint"),
1645 fHM->
H1(
"fhPhotonEPmtHit")},
1646 {
"Sens plane Z+",
"Sens plane Z-",
"PMT point",
"PMT hit"});
1649 fHM->
H1(
"fhLambdaPlaneZ-"),
1650 fHM->
H1(
"fhLambdaPmtPoint"),
1651 fHM->
H1(
"fhLambdaPmtHit")},
1652 {
"Sens plane Z+",
"Sens plane Z-",
"PMT Point",
"PMT hit"});
1659 "fhPointsXYZ;X [cm];Y [cm];Z [cm];Yield",
1670 "fhHitsXYZ;X [cm];Y [cm];Z [cm];Yield",
1680 vector<Int_t> pixels =
1686 fHM->
CreateCanvas(
"richgeo_pixels_xy",
"richgeo_pixels_xy", 1500, 1500);
1687 c->SetGrid(
true,
true);
1688 TH2D* pad =
new TH2D(
1689 "richgeo_pixels_xy",
";X [cm];Y [cm]", 1, -120, 120, 1, -210, 210);
1691 pad->SetStats(
false);
1698 fHM->
CreateCanvas(
"richgeo_pixels_xz",
"richgeo_pixels_xz", 1500, 1500);
1699 c->SetGrid(
true,
true);
1700 TH2D* pad =
new TH2D(
1701 "richgeo_pixels_xz",
";Z [cm];X [cm]", 1, 200, 250, 1, -120., 120.);
1702 pad->SetStats(
false);
1708 fHM->
CreateCanvas(
"richgeo_pixels_yz",
"richgeo_pixels_yz", 1500, 1500);
1709 TH2D* pad =
new TH2D(
1710 "richgeo_pixels_yz",
";Z [cm];Y [cm]", 1, 200, 250, 1, -220, 220);
1711 pad->SetStats(
false);
1719 c->SetGrid(
true,
true);
1721 new TH2D(
"richgeo_pmts_xy",
";X [cm];Y [cm]", 1, -120, 120, 1, -210, 210);
1722 pad->SetStats(
false);
1728 for (
unsigned int i = 0;
i < pixels.size();
i++) {
1731 TVector3 inPos(pixelData->
fX, pixelData->
fY, pixelData->
fZ);
1735 fHM->
H3(
"fhPointsXYZ")->Fill(inPos.X(), inPos.Y(), inPos.Z());
1736 fHM->
H3(
"fhHitsXYZ")->Fill(outPos.X(), outPos.Y(), outPos.Z());
1741 "richgeo_pixels_points_xyz",
"richgeo_pixels_points_xyz", 1500, 1500);
1742 fHM->
H3(
"fhPointsXYZ")->Draw();
1747 "richgeo_pixels_hits_xyz",
"richgeo_pixels_hits_xyz", 1500, 1500);
1748 fHM->
H3(
"fhHitsXYZ")->Draw();
1753 const vector<Int_t>& ids,
1754 Bool_t isDrawPixel) {
1755 for (
unsigned int i = 0;
i < ids.size();
i++) {
1757 Double_t boxHalfSize = 0.0;
1761 inPos.SetXYZ(pixelData->
fX, pixelData->
fY, pixelData->
fZ);
1766 inPos.SetXYZ(pmtData->
fX, pmtData->
fY, pmtData->
fZ);
1767 boxHalfSize = 0.5 * pmtData->
fHeight;
1771 TBox* boxOut =
nullptr;
1773 if (coordinates ==
string(
"xy"))
1774 boxOut =
new TBox(outPos.X() - boxHalfSize,
1775 outPos.Y() - boxHalfSize,
1776 outPos.X() + boxHalfSize,
1777 outPos.Y() + boxHalfSize);
1778 if (coordinates ==
string(
"zx"))
1779 boxOut =
new TBox(outPos.Z() - boxHalfSize,
1780 outPos.X() - boxHalfSize,
1781 outPos.Z() + boxHalfSize,
1782 outPos.X() + boxHalfSize);
1783 if (coordinates ==
string(
"zy"))
1784 boxOut =
new TBox(outPos.Z() - boxHalfSize,
1785 outPos.Y() - boxHalfSize,
1786 outPos.Z() + boxHalfSize,
1787 outPos.Y() + boxHalfSize);
1788 if (boxOut != NULL) {
1790 boxOut->SetFillColor(kBlue);
1792 boxOut->SetFillStyle(0);
1793 boxOut->SetLineColor(kBlue);
1798 TBox* boxIn =
nullptr;
1799 if (coordinates ==
string(
"xy"))
1800 boxIn =
new TBox(inPos.X() - boxHalfSize,
1801 inPos.Y() - boxHalfSize,
1802 inPos.X() + boxHalfSize,
1803 inPos.Y() + boxHalfSize);
1804 if (coordinates ==
string(
"zx"))
1805 boxIn =
new TBox(inPos.Z() - boxHalfSize,
1806 inPos.X() - boxHalfSize,
1807 inPos.Z() + boxHalfSize,
1808 inPos.X() + boxHalfSize);
1809 if (coordinates ==
string(
"zy"))
1810 boxIn =
new TBox(inPos.Z() - boxHalfSize,
1811 inPos.Y() - boxHalfSize,
1812 inPos.Z() + boxHalfSize,
1813 inPos.Y() + boxHalfSize);
1814 if (boxIn != NULL) {
1816 boxIn->SetFillColor(kRed);
1818 boxIn->SetFillStyle(0);
1819 boxIn->SetLineColor(kRed);
1831 TDirectory* oldir = gDirectory;
1832 TFile* outFile = FairRootManager::Instance()->GetOutFile();
1833 if (outFile != NULL) {
1837 gDirectory->cd(oldir->GetPath());
1842 if (histAcc->GetEntries() == 0) {
1846 100. * Double_t(histRec->GetEntries()) / Double_t(histAcc->GetEntries());
1852 const string& outputDir) {
1855 if (
fHM != NULL)
delete fHM;
1858 TFile* file =
new TFile(fileName.c_str());