CbmRichUrqmdTest.cxx 24.17 KiB
/* Copyright (C) 2012-2024 UGiessen/JINR-LIT, Giessen/Dubna
SPDX-License-Identifier: GPL-3.0-only
Authors: Semen Lebedev [committer], Andrey Lebedev, Martin Beyer */
#include "CbmRichUrqmdTest.h"
#include "CbmDigiManager.h"
#include "CbmDrawHist.h"
#include "CbmHistManager.h"
#include "CbmMCDataArray.h"
#include "CbmMCDataManager.h"
#include "CbmMCEventList.h"
#include "CbmMCTrack.h"
#include "CbmMatchRecoToMC.h"
#include "CbmRichDetectorData.h"
#include "CbmRichDigi.h"
#include "CbmRichDigiMapManager.h"
#include "CbmRichDraw.h"
#include "CbmRichGeoManager.h"
#include "CbmRichHit.h"
#include "CbmRichPoint.h"
#include "CbmRichRing.h"
#include "CbmRichUtil.h"
#include "CbmTrackMatchNew.h"
#include "CbmUtils.h"
#include "FairMCPoint.h"
#include "FairTrackParam.h"
#include "TCanvas.h"
#include "TClonesArray.h"
#include "TH1.h"
#include "TH1D.h"
#include "TStyle.h"
#include <TFile.h>
#include <iostream>
#include <sstream>
#include <string>
using namespace std;
using namespace Cbm;
CbmRichUrqmdTest::CbmRichUrqmdTest() : FairTask("CbmRichUrqmdTest") {}
CbmRichUrqmdTest::~CbmRichUrqmdTest() {}
InitStatus CbmRichUrqmdTest::Init()
{
fMcTracks = InitOrFatalMc("MCTrack", "CbmRichUrqmdTest::Init");
fRichPoints = InitOrFatalMc("RichPoint", "CbmRichUrqmdTest::Init");
fRichHits = GetOrFatal<TClonesArray>("RichHit", "CbmRichUrqmdTest::Init");
fRichRings = GetOrFatal<TClonesArray>("RichRing", "CbmRichUrqmdTest::Init");
fRichRingMatches = GetOrFatal<TClonesArray>("RichRingMatch", "CbmRichUrqmdTest::Init");
fRichProjections = GetOrFatal<TClonesArray>("RichProjection", "CbmRichUrqmdTest::Init");
fEventList = GetOrFatal<CbmMCEventList>("MCEventList.", "CbmRichUrqmdTest::Init");
fDigiMan = CbmDigiManager::Instance();
fDigiMan->Init();
fVertexZStsSlices = {make_pair(0., 5.), make_pair(5., 15.), make_pair(15., 25.), make_pair(25., 35.),
make_pair(35., 45.), make_pair(45., 55.), make_pair(55., 65.), make_pair(65., 75.),
make_pair(75., 85.), make_pair(85., 95.), make_pair(95., 105.)};
InitHistograms();
return kSUCCESS;
}
void CbmRichUrqmdTest::Exec(Option_t* /*option*/){
fEventNum++;
LOG(info) << "CbmRichUrqmdTest, event No. " << fEventNum;
FillRichRingNofHits();
NofRings();
NofHitsAndPoints();
NofProjections();
Vertex();
PmtXYSource();
}
void CbmRichUrqmdTest::FillRichRingNofHits()
{
fNofHitsInRingMap.clear();
int nofRichHits = fRichHits->GetEntriesFast();
for (int iHit = 0; iHit < nofRichHits; iHit++) {
const CbmRichHit* hit = static_cast<CbmRichHit*>(fRichHits->At(iHit));
if (nullptr == hit) continue;
vector<CbmLink> motherIds = CbmMatchRecoToMC::GetMcTrackMotherIdsForRichHit(fDigiMan, hit, fRichPoints, fMcTracks);
for (const auto& motherId : motherIds) {
fNofHitsInRingMap[motherId]++;
}
}
}
void CbmRichUrqmdTest::InitHistograms()
{
fHM = new CbmHistManager();
fHM->Create1<TH1D>("fh_vertex_z", "fh_vertex_z;Z [cm];# vertices/ev.", 400, -50., 350);
fHM->Create1<TH1D>("fh_vertex_z_sts", "fh_vertex_z_sts;Z [cm];# vertices/ev.", 320, -50., 110.);
fHM->Create2<TH2D>("fh_vertex_xy", "fh_vertex_xy;X [cm];Y [cm];# vertices/ev.", 100, -200., 200., 100, -200., 200.);
fHM->Create2<TH2D>("fh_vertex_zy", "fh_vertex_zy;Z [cm];Y [cm];# vertices/ev.", 400, -50., 350, 100, -200., 200.);
fHM->Create2<TH2D>("fh_vertex_zx", "fh_vertex_zx;Z [cm];X [cm];# vertices/ev.", 400, -50., 350, 100, -200., 200.);
vector<string> vertexZTypes {"z60_140", "z140_330", "z140_190"};
for (const string& t : vertexZTypes) {
fHM->Create2<TH2D>("fh_vertex_xy_" + t, "fh_vertex_xy_" + t + ";X [cm];Y [cm];# vertices/ev.", 100, -200., 200.,
100, -200., 200.);
}
for (auto pair : fVertexZStsSlices) {
string name = "fh_vertex_xy_z" + to_string(pair.first) + "_" + to_string(pair.second);
fHM->Create2<TH2D>(name, name + ";x [cm];y [cm];# vertices/ev.", 100, -100., 100., 100, -100., 100.);
}
vector<string> nofRingsTypes {"1hit", "7hits", "prim_1hit", "prim_7hits", "target_1hit", "target_7hits"};
for (const string& t : nofRingsTypes) {
double nofBins = (t == "1hit" || t == "7hits") ? 30 : 100;
fHM->Create1<TH1D>("fh_nof_rings_" + t, "fh_nof_rings_" + t + ";# particles/ev.;Yield", nofBins, -.5,
nofBins - 0.5);
}
vector<string> momP {"fh_secel_mom", "fh_gamma_target_mom", "fh_gamma_nontarget_mom",
"fh_pi_mom", "fh_kaon_mom", "fh_mu_mom"};
for (const string& t : momP) {
fHM->Create1<TH1D>(t, t + ";P [GeV/c];Number per event", 100, 0., 20);
}
fHM->Create1<TH1D>("fh_nof_points_per_event_src", "fh_nof_points_per_event_src;Particle;# MC points/ev.", 7, .5, 7.5);
fHM->Create1<TH1D>("fh_nof_hits_per_event", "fh_nof_hits_per_event;# hits/ev.;Yield", 200, 0, 2000);
fHM->Create1<TH1D>("fh_nof_points_per_event", "fh_nof_points_per_event;# points/ev.;Yield", 200, 0, 10000);
fHM->Create1<TH1D>("fh_nof_hits_per_pmt", "fh_nof_hits_per_pmt;# hits/PMT;% of total", 65, -0.5, 64.5);
vector<double> xPmtBins = CbmRichUtil::GetPmtHistXbins();
vector<double> yPmtBins = CbmRichUtil::GetPmtHistYbins();
// before drawing this histogram must be normalized by 1/64
fHM->Create2<TH2D>("fh_hitrate_xy", "fh_hitrate_xy;X [cm];Y [cm];# hits/pixel/s", xPmtBins, yPmtBins);
fHM->Create2<TH2D>("fh_hits_xy", "fh_hits_xy;X [cm];Y [cm];# hits/PMT/ev.", xPmtBins, yPmtBins);
vector<string> pointXYTypes {"", "_pions", "_gamma_target", "_gamma_nontarget"};
for (const string& t : pointXYTypes) {
fHM->Create2<TH2D>("fh_points_xy" + t, "fh_points_xy" + t + ";X [cm];Y [cm];# MC points/PMT/ev.", xPmtBins,
yPmtBins);
}
vector<string> skippedPmtTypes {"10", "20", "30"};
for (const string& t : skippedPmtTypes) {
fHM->Create2<TH2D>("fh_skipped_pmt_xy_" + t,
"fh_skipped_pmt_xy_" + t + ";X [cm];Y [cm];# skipped PMTs (>" + t + " hits) [%]", xPmtBins,
yPmtBins);
}
fHM->Create1<TH1D>("fh_nof_proj_per_event", "fh_nof_proj_per_event;# tracks/ev.;Yield", 600, -.5, 599.5);
fHM->Create2<TH2D>("fh_proj_xy", "fh_proj_xy;X [cm];Y [cm];# tracks/cm^{2}/ev.", 240, -120, 120, 420, -210, 210);
}
void CbmRichUrqmdTest::NofRings()
{
int nofRings = fRichRings->GetEntriesFast();
int nRings1hit = 0, nRings7hits = 0;
int nRingsPrim1hit = 0, nRingsPrim7hits = 0;
int nRingsTarget1hit = 0, nRingsTarget7hits = 0;
for (int iR = 0; iR < nofRings; iR++) {
const CbmRichRing* ring = static_cast<CbmRichRing*>(fRichRings->At(iR));
if (ring == nullptr) continue;
const CbmTrackMatchNew* ringMatch = static_cast<CbmTrackMatchNew*>(fRichRingMatches->At(iR));
if (ringMatch == nullptr || ringMatch->GetNofLinks() < 1) continue;
auto matchedLink = ringMatch->GetMatchedLink();
const CbmMCTrack* mcTrack = static_cast<CbmMCTrack*>(fMcTracks->Get(matchedLink));
if (mcTrack == nullptr) continue;
int motherId = mcTrack->GetMotherId();
int pdgAbs = std::abs(mcTrack->GetPdgCode());
double mom = mcTrack->GetP();
TVector3 vert;
mcTrack->GetStartVertex(vert);
double dZ = vert.Z();
if (motherId == -1 && pdgAbs == 11) continue; // do not calculate embedded electrons
int nofHits = ring->GetNofHits();
if (nofHits >= 1) nRings1hit++;
if (nofHits >= fMinNofHits) nRings7hits++;
if (motherId == -1 && nofHits >= 1) nRingsPrim1hit++;
if (motherId == -1 && nofHits >= fMinNofHits) nRingsPrim7hits++;
if (dZ < 0.1 && nofHits >= 1) nRingsTarget1hit++;
if (dZ < 0.1 && nofHits >= fMinNofHits) nRingsTarget7hits++;
if (nofHits >= 1) {
if (motherId != -1) {
int motherPdg = -1;
const CbmMCTrack* mother = static_cast<CbmMCTrack*>(fMcTracks->Get(matchedLink));
if (mother != nullptr) motherPdg = mother->GetPdgCode();
if (motherId != -1 && pdgAbs == 11 && motherPdg != 22) fHM->H1("fh_secel_mom")->Fill(mom);
if (motherId != -1 && pdgAbs == 11 && motherPdg == 22) {
if (dZ < 0.1) { fHM->H1("fh_gamma_target_mom")->Fill(mom); }
else {
fHM->H1("fh_gamma_nontarget_mom")->Fill(mom);
}
}
} if (pdgAbs == 211) fHM->H1("fh_pi_mom")->Fill(mom);
if (pdgAbs == 321) fHM->H1("fh_kaon_mom")->Fill(mom);
if (pdgAbs == 13) fHM->H1("fh_mu_mom")->Fill(mom);
}
}
fHM->H1("fh_nof_rings_1hit")->Fill(nRings1hit);
fHM->H1("fh_nof_rings_7hits")->Fill(nRings7hits);
fHM->H1("fh_nof_rings_prim_1hit")->Fill(nRingsPrim1hit);
fHM->H1("fh_nof_rings_prim_7hits")->Fill(nRingsPrim7hits);
fHM->H1("fh_nof_rings_target_1hit")->Fill(nRingsTarget1hit);
fHM->H1("fh_nof_rings_target_7hits")->Fill(nRingsTarget7hits);
}
void CbmRichUrqmdTest::NofHitsAndPoints()
{
int nofHits = fRichHits->GetEntriesFast();
fHM->H1("fh_nof_hits_per_event")->Fill(nofHits);
map<int, int> digisPerPmtMap;
for (int iH = 0; iH < nofHits; iH++) {
const CbmRichHit* hit = static_cast<CbmRichHit*>(fRichHits->At(iH));
if (hit == nullptr) continue;
fHM->H2("fh_hits_xy")->Fill(hit->GetX(), hit->GetY());
fHM->H2("fh_hitrate_xy")->Fill(hit->GetX(), hit->GetY());
int digiInd = hit->GetRefId();
const CbmRichDigi* richDigi = fDigiMan->Get<CbmRichDigi>(digiInd);
if (richDigi == nullptr) continue;
CbmRichPixelData* pixelData = CbmRichDigiMapManager::GetInstance().GetPixelDataByAddress(richDigi->GetAddress());
if (pixelData == nullptr) continue;
int pmtId = pixelData->fPmtId;
digisPerPmtMap[pmtId]++;
}
for (auto const& it : digisPerPmtMap) {
int pmtId = it.first;
int nofDigis = it.second;
CbmRichPmtData* pmtData = CbmRichDigiMapManager::GetInstance().GetPmtDataById(pmtId);
TVector3 inPos(pmtData->fX, pmtData->fY, pmtData->fZ);
TVector3 outPos;
CbmRichGeoManager::GetInstance().RotatePoint(&inPos, &outPos);
if (nofDigis > 10) fHM->H2("fh_skipped_pmt_xy_10")->Fill(outPos.X(), outPos.Y());
if (nofDigis > 20) fHM->H2("fh_skipped_pmt_xy_20")->Fill(outPos.X(), outPos.Y());
if (nofDigis > 30) fHM->H2("fh_skipped_pmt_xy_30")->Fill(outPos.X(), outPos.Y());
}
vector<int> allPmtIds = CbmRichDigiMapManager::GetInstance().GetPmtIds();
for (int pmtId : allPmtIds) {
fHM->H1("fh_nof_hits_per_pmt")->Fill(digisPerPmtMap[pmtId]);
}
int nofEvents = fEventList->GetNofEvents();
for (int iE = 0; iE < nofEvents; iE++) {
int fileId = fEventList->GetFileIdByIndex(iE);
int eventId = fEventList->GetEventIdByIndex(iE);
int nofPoints = fRichPoints->Size(fileId, eventId);
fHM->H1("fh_nof_points_per_event")->Fill(nofPoints);
for (int i = 0; i < nofPoints; i++) {
const CbmRichPoint* point = static_cast<CbmRichPoint*>(fRichPoints->Get(fileId, eventId, i));
if (point == nullptr) continue;
fHM->H1("fh_nof_points_per_event_src")->Fill(1);
int mcPhotonTrackId = point->GetTrackID();
if (mcPhotonTrackId < 0) continue;
const CbmMCTrack* mcPhotonTrack = static_cast<CbmMCTrack*>(fMcTracks->Get(fileId, eventId, mcPhotonTrackId));
if (mcPhotonTrack == nullptr) continue;
int motherPhotonId = mcPhotonTrack->GetMotherId();
if (motherPhotonId < 0) continue;
const CbmMCTrack* mcTrack = static_cast<CbmMCTrack*>(fMcTracks->Get(fileId, eventId, motherPhotonId)); if (mcTrack == nullptr) continue;
int motherId = mcTrack->GetMotherId();
int pdgAbs = std::abs(mcTrack->GetPdgCode());
TVector3 vert;
mcTrack->GetStartVertex(vert);
double dZ = vert.Z();
if (motherId == -1 && pdgAbs == 11) continue; // do not calculate embedded electrons
if (motherId != -1) {
int motherPdg = -1;
const CbmMCTrack* mother = static_cast<CbmMCTrack*>(fMcTracks->Get(fileId, eventId, motherId));
if (mother != nullptr) motherPdg = mother->GetPdgCode();
if (motherId != -1 && pdgAbs == 11 && motherPdg != 22) fHM->H1("fh_nof_points_per_event_src")->Fill(2);
if (motherId != -1 && pdgAbs == 11 && motherPdg == 22) {
if (dZ < 0.1) { fHM->H1("fh_nof_points_per_event_src")->Fill(3); }
else {
fHM->H1("fh_nof_points_per_event_src")->Fill(4);
}
}
}
if (pdgAbs == 211) fHM->H1("fh_nof_points_per_event_src")->Fill(5);
if (pdgAbs == 321) fHM->H1("fh_nof_points_per_event_src")->Fill(6);
if (pdgAbs == 13) fHM->H1("fh_nof_points_per_event_src")->Fill(7);
}
}
}
void CbmRichUrqmdTest::PmtXYSource()
{
int nofEvents = fEventList->GetNofEvents();
for (int iE = 0; iE < nofEvents; iE++) {
int fileId = fEventList->GetFileIdByIndex(iE);
int eventId = fEventList->GetEventIdByIndex(iE);
int nofPoints = fRichPoints->Size(fileId, eventId);
for (int i = 0; i < nofPoints; i++) {
const CbmRichPoint* point = static_cast<CbmRichPoint*>(fRichPoints->Get(fileId, eventId, i));
if (point == nullptr) continue;
int iMCTrack = point->GetTrackID();
const CbmMCTrack* track = static_cast<CbmMCTrack*>(fMcTracks->Get(fileId, eventId, iMCTrack));
if (track == nullptr) continue;
int iMother = track->GetMotherId();
if (iMother == -1) continue;
const CbmMCTrack* track2 = static_cast<CbmMCTrack*>(fMcTracks->Get(fileId, eventId, iMother));
if (track2 == nullptr) continue;
int pdgAbs = std::abs(track2->GetPdgCode());
int motherId = track2->GetMotherId();
TVector3 inPos(point->GetX(), point->GetY(), point->GetZ());
TVector3 outPos;
CbmRichGeoManager::GetInstance().RotatePoint(&inPos, &outPos);
fHM->H2("fh_points_xy")->Fill(outPos.X(), outPos.Y());
if (motherId != -1) {
int motherPdg = -1;
const CbmMCTrack* mother = static_cast<CbmMCTrack*>(fMcTracks->Get(fileId, eventId, motherId));
if (mother != nullptr) motherPdg = mother->GetPdgCode();
TVector3 vert;
track2->GetStartVertex(vert);
if (motherId != -1 && pdgAbs == 11 && motherPdg == 22) {
if (vert.Z() < 0.1) { fHM->H2("fh_points_xy_gamma_target")->Fill(outPos.X(), outPos.Y()); }
else {
fHM->H2("fh_points_xy_gamma_nontarget")->Fill(outPos.X(), outPos.Y());
}
}
} if (pdgAbs == 211) fHM->H2("fh_points_xy_pions")->Fill(outPos.X(), outPos.Y());
}
}
}
void CbmRichUrqmdTest::NofProjections()
{
if (fRichProjections == nullptr) return;
int nofProj = fRichProjections->GetEntriesFast();
int nofGoodProj = 0;
for (int i = 0; i < nofProj; i++) {
FairTrackParam* proj = (FairTrackParam*) fRichProjections->At(i);
if (proj == nullptr) continue;
fHM->H2("fh_proj_xy")->Fill(proj->GetX(), proj->GetY());
if (proj->GetX() != 0 && proj->GetY() != 0) nofGoodProj++;
}
fHM->H1("fh_nof_proj_per_event")->Fill(nofGoodProj);
}
void CbmRichUrqmdTest::Vertex()
{
int nofEvents = fEventList->GetNofEvents();
for (int iE = 0; iE < nofEvents; iE++) {
int fileId = fEventList->GetFileIdByIndex(iE);
int eventId = fEventList->GetEventIdByIndex(iE);
int nMcTracks = fMcTracks->Size(fileId, eventId);
for (int i = 0; i < nMcTracks; i++) {
//At least one hit in RICH
CbmLink val(1., i, eventId, fileId);
if (fNofHitsInRingMap[val] < 1) continue;
const CbmMCTrack* mctrack = static_cast<CbmMCTrack*>(fMcTracks->Get(val));
TVector3 v;
mctrack->GetStartVertex(v);
fHM->H1("fh_vertex_z")->Fill(v.Z());
fHM->H1("fh_vertex_z_sts")->Fill(v.Z());
fHM->H2("fh_vertex_xy")->Fill(v.X(), v.Y());
fHM->H2("fh_vertex_zy")->Fill(v.Z(), v.Y());
fHM->H2("fh_vertex_zx")->Fill(v.Z(), v.X());
if (v.Z() >= 60 && v.Z() <= 140) fHM->H2("fh_vertex_xy_z60_140")->Fill(v.X(), v.Y());
if (v.Z() >= 140 && v.Z() <= 330) fHM->H2("fh_vertex_xy_z140_330")->Fill(v.X(), v.Y());
if (v.Z() >= 140 && v.Z() <= 190) fHM->H2("fh_vertex_xy_z140_190")->Fill(v.X(), v.Y());
for (auto pair : fVertexZStsSlices) {
string name = "fh_vertex_xy_z" + to_string(pair.first) + "_" + to_string(pair.second);
if (v.Z() > pair.first && v.Z() <= pair.second) { fHM->H2(name)->Fill(v.X(), v.Y()); }
}
}
}
}
void CbmRichUrqmdTest::DrawHist()
{
cout.precision(4);
SetDefaultDrawStyle();
{
fHM->Scale("fh_vertex_z", 1. / fEventNum);
fHM->CreateCanvas("richurqmd_vertex_z", "richurqmd_vertex_z", 1000, 1000);
DrawH1(fHM->H1("fh_vertex_z"), kLinear, kLog, "hist");
}
{
fHM->Scale("fh_vertex_z_sts", 1. / fEventNum);
fHM->CreateCanvas("richurqmd_vertex_z_sts", "richurqmd_vertex_z_sts", 1000, 1000);
DrawH1(fHM->H1("fh_vertex_z_sts"), kLinear, kLog, "hist");
}
{
fHM->Scale("fh_vertex_xy", 1. / fEventNum);
fHM->Scale("fh_vertex_zy", 1. / fEventNum);
fHM->Scale("fh_vertex_zx", 1. / fEventNum);
TCanvas* c = fHM->CreateCanvas("richurqmd_vertex_xyz", "richurqmd_vertex_xyz", 1800, 600);
c->Divide(3, 1);
c->cd(1);
DrawH2(fHM->H2("fh_vertex_xy"), kLinear, kLinear, kLog);
c->cd(2);
DrawH2(fHM->H2("fh_vertex_zy"), kLinear, kLinear, kLog);
c->cd(3);
DrawH2(fHM->H2("fh_vertex_zx"), kLinear, kLinear, kLog);
}
{
gStyle->SetOptTitle(1);
TCanvas* c = fHM->CreateCanvas("richurqmd_vertex_sts_xyz", "richurqmd_vertex_sts_xyz", 1600, 1200);
c->Divide(4, 3);
int i = 1;
for (auto pair : fVertexZStsSlices) {
string name = "fh_vertex_xy_z" + to_string(pair.first) + "_" + to_string(pair.second);
fHM->Scale(name, 1. / fEventNum);
c->cd(i++);
DrawH2(fHM->H2(name), kLinear, kLinear, kLog);
DrawTextOnPad(to_string(pair.first) + " cm < Z < " + to_string(pair.second) + " cm", 0.3, 0.9, 0.7, 0.98);
}
gStyle->SetOptTitle(0);
}
vector<string> vertexZTypes {"z60_140", "z140_330", "z140_190"};
for (const string& t : vertexZTypes) {
string name = "richurqmd_vertex_xy_" + t;
fHM->Scale("fh_vertex_xy_" + t, 1. / fEventNum);
fHM->CreateCanvas(name, name, 1000, 1000);
DrawH2(fHM->H2("fh_vertex_xy_" + t));
}
{
fHM->Scale("fh_nof_points_per_event_src", 1. / fEventNum);
fHM->CreateCanvas("richurqmd_nof_points_per_event_src", "richurqmd_nof_points_per_event_src", 1000, 1000);
//gStyle->SetPaintTextFormat("4.1f");
string labels[7] = {"all",
"e^{#pm}_{sec} other",
"e^{#pm}_{target} from #gamma",
"e^{#pm}_{not target} from #gamma",
"#pi^{#pm}",
"K^{#pm}",
"#mu^{#pm}"};
for (int i = 1; i <= 7; i++) {
fHM->H1("fh_nof_points_per_event_src")->GetXaxis()->SetBinLabel(i, labels[i - 1].c_str());
}
fHM->H1("fh_nof_points_per_event_src")->GetXaxis()->SetLabelSize(0.05);
//fHM->H1("fh_nof_points_per_event_src")->SetMarkerSize(0.15);
DrawH1(fHM->H1("fh_nof_points_per_event_src"), kLinear, kLog, "hist text0");
}
{
vector<string> nofRingsTypes {"", "_prim", "_target"};
for (const string& t : nofRingsTypes) {
string cName = "richurqmd_nof_rings" + t;
string h1Name = "fh_nof_rings" + t + "_1hit";
string h7Name = "fh_nof_rings" + t + "_7hits";
fHM->CreateCanvas(cName, cName, 1000, 1000);
fHM->NormalizeToIntegral(h1Name);
fHM->NormalizeToIntegral(h7Name);
stringstream ss1, ss2;
ss1 << "At least 1 hit (" << fHM->H1(h1Name)->GetMean() << ")";
ss2 << "At least 7 hits (" << fHM->H1(h7Name)->GetMean() << ")"; DrawH1({fHM->H1(h1Name), fHM->H1(h7Name)}, {ss1.str(), ss2.str()}, kLinear, kLinear, true, 0.4, 0.85, 0.99, 0.99,
"hist");
}
}
{
fHM->CreateCanvas("richurqmd_sources_mom", "richurqmd_sources_mom", 1000, 1000);
fHM->Scale("fh_gamma_target_mom", 1. / fEventNum);
fHM->Scale("fh_gamma_nontarget_mom", 1. / fEventNum);
fHM->Scale("fh_secel_mom", 1. / fEventNum);
fHM->Scale("fh_pi_mom", 1. / fEventNum);
fHM->Scale("fh_kaon_mom", 1. / fEventNum);
fHM->Scale("fh_mu_mom", 1. / fEventNum);
stringstream ss1, ss2, ss3, ss4, ss5, ss6;
ss1 << "e^{#pm}_{target} from #gamma (" << fHM->H1("fh_gamma_target_mom")->GetEntries() / fEventNum << ")";
ss2 << "e^{#pm}_{not target} from #gamma (" << fHM->H1("fh_gamma_nontarget_mom")->GetEntries() / fEventNum << ")";
ss3 << "e^{#pm}_{sec} other (" << fHM->H1("fh_secel_mom")->GetEntries() / fEventNum << ")";
ss4 << "#pi^{#pm} (" << fHM->H1("fh_pi_mom")->GetEntries() / fEventNum << ")";
ss5 << "K^{#pm} (" << fHM->H1("fh_kaon_mom")->GetEntries() / fEventNum << ")";
ss6 << "#mu^{#pm} (" << fHM->H1("fh_mu_mom")->GetEntries() / fEventNum << ")";
DrawH1(fHM->H1Vector({"fh_gamma_target_mom", "fh_gamma_nontarget_mom", "fh_secel_mom", "fh_pi_mom", "fh_kaon_mom",
"fh_mu_mom"}),
{ss1.str(), ss2.str(), ss3.str(), ss4.str(), ss5.str(), ss6.str()}, kLinear, kLog, true, 0.5, 0.7, 0.99,
0.99, "hist");
}
{
TCanvas* c = fHM->CreateCanvas("richurqmd_hits_xy", "richurqmd_hits_xy", 1000, 1000);
TH2* clone = fHM->H2Clone("fh_hits_xy");
clone->Scale(1. / (fEventNum));
CbmRichDraw::DrawPmtH2(clone, c, true);
}
{
TCanvas* c = fHM->CreateCanvas("richurqmd_occupancy_xy", "richurqmd_occupancy_xy", 1000, 1000);
TH2* clone = fHM->H2Clone("fh_hits_xy");
clone->GetZaxis()->SetTitle("Occupancy:# hits/PMT/ev./64 [%]");
clone->Scale(100. / (fEventNum * 64.));
CbmRichDraw::DrawPmtH2(clone, c, true);
}
{
vector<string> skippedPmtTypes {"10", "20", "30"};
for (const string& t : skippedPmtTypes) {
string name = "richurqmd_skipped_pmt_xy_" + t;
TCanvas* c = fHM->CreateCanvas(name, name, 1000, 1000);
fHM->Scale("fh_skipped_pmt_xy_" + t, 100. / fEventNum);
CbmRichDraw::DrawPmtH2(fHM->H2("fh_skipped_pmt_xy_" + t), c, true);
}
}
{
fHM->CreateCanvas("richurqmd_nof_hits_per_pmt", "richurqmd_nof_hits_per_pmt", 1000, 1000);
fHM->NormalizeToIntegral("fh_nof_hits_per_pmt");
DrawH1(fHM->H1("fh_nof_hits_per_pmt"), kLinear, kLog, "hist");
fHM->H1("fh_nof_hits_per_pmt")->SetStats(true);
}
{
vector<string> pointXYTypes {"", "_pions", "_gamma_target", "_gamma_nontarget"};
for (const string& t : pointXYTypes) {
string name = "richurqmd_points_xy" + t;
TCanvas* c = fHM->CreateCanvas(name, name, 1000, 1000);
fHM->Scale("fh_points_xy" + t, 1. / fEventNum);
CbmRichDraw::DrawPmtH2(fHM->H2("fh_points_xy" + t), c, true);
}
}
{
fHM->CreateCanvas("richurqmd_nof_hits_per_event", "richurqmd_nof_hits_per_event", 1000, 1000); fHM->NormalizeToIntegral("fh_nof_hits_per_event");
DrawH1(fHM->H1("fh_nof_hits_per_event"), kLinear, kLinear, "hist");
fHM->H1("fh_nof_hits_per_event")->SetStats(true);
cout << "Mean number of hits per event = " << fHM->H1("fh_nof_hits_per_event")->GetMean() << endl;
}
{
fHM->CreateCanvas("richurqmd_nof_points_per_event", "richurqmd_nof_points_per_event", 1000, 1000);
fHM->NormalizeToIntegral("fh_nof_points_per_event");
DrawH1(fHM->H1("fh_nof_points_per_event"), kLinear, kLinear, "hist");
fHM->H1("fh_nof_points_per_event")->SetStats(true);
cout << "Mean number of points per event = " << fHM->H1("fh_nof_points_per_event")->GetMean() << endl;
}
{
TCanvas* c = fHM->CreateCanvas("richurqmd_hitrate_xy", "richurqmd_hitrate_xy", 1000, 1000);
fHM->Scale("fh_hitrate_xy", 1e7 / (fEventNum * 64.));
CbmRichDraw::DrawPmtH2(fHM->H2("fh_hitrate_xy"), c, true);
//DrawH2(fHM->H2("fh_hitrate_xy"));
}
{
TCanvas* c = fHM->CreateCanvas("richurqmd_proj_xy", "richurqmd_proj_xy", 1000, 1000);
fHM->Scale("fh_proj_xy", 1. / fEventNum);
CbmRichDraw::DrawPmtH2(fHM->H2("fh_proj_xy"), c);
}
{
fHM->CreateCanvas("richurqmd_nof_proj_per_event", "richurqmd_nof_proj_per_event", 1000, 1000);
fHM->Scale("fh_nof_proj_per_event", 1. / fEventNum);
DrawH1(fHM->H1("fh_nof_proj_per_event"), kLinear, kLinear, "hist");
fHM->H1("fh_nof_proj_per_event")->SetStats(true);
cout << "Number of track projections per event = " << fHM->H1("fh_nof_proj_per_event")->GetMean() << endl;
}
}
void CbmRichUrqmdTest::Finish()
{
TDirectory* oldir = gDirectory;
TFile* outFile = FairRootManager::Instance()->GetOutFile();
if (outFile != nullptr) {
outFile->mkdir(GetName());
outFile->cd(GetName());
fHM->WriteToFile();
}
DrawHist();
fHM->SaveCanvasToImage(fOutputDir);
fHM->Clear();
gDirectory->cd(oldir->GetPath());
}
ClassImp(CbmRichUrqmdTest)