Commit e7f1e0ef authored by Semen Lebedev's avatar Semen Lebedev
Browse files

Updates in RichUrqmdTestQa

parent f933e3a8
......@@ -88,6 +88,21 @@ public:
Add(name, h);
}
/**
* \brief Helper function for creation of 1-dimensional histograms and profiles.
* Template argument is a real object type that has to be created, for example,
* Create1<TH1F>("name", "title", binsX);
* \param[in] name Object name.
* \param[in] title Object title.
* \param[in] binsX array of low-edges for each bin in X
*/
template<class T>
void Create1(const std::string& name, const std::string& title, const std::vector<Double_t>& binsX)
{
T* h = new T(name.c_str(), title.c_str(), binsX.size() - 1, binsX.data());
Add(name, h);
}
/**
* \brief Helper function for creation of 2-dimensional histograms and profiles.
* Template argument is a real object type that has to be created, for example,
......@@ -109,6 +124,23 @@ public:
Add(name, h);
}
/**
* \brief Helper function for creation of 2-dimensional histograms and profiles.
* Template argument is a real object type that has to be created, for example,
* Create2<TH2F>("name", "title", binsX, binsY);
* \param[in] name Object name.
* \param[in] title Object title.
* \param[in] binsX array of low-edges for each bin in X
* \param[in] binsY array of low-edges for each bin in Y
*/
template<class T>
void Create2(const std::string& name, const std::string& title, const std::vector<Double_t>& binsX,
const std::vector<Double_t>& binsY)
{
T* h = new T(name.c_str(), title.c_str(), binsX.size() - 1, binsX.data(), binsY.size() - 1, binsY.data());
Add(name, h);
}
/**
* \brief Helper function for creation of 3-dimensional histograms and profiles.
* Template argument is a real object type that has to be created, for example,
......
......@@ -25,6 +25,14 @@ using std::vector;
namespace Cbm
{
CbmMCDataArray* InitOrFatalMc(const std::string& objName, const std::string& description)
{
CbmMCDataManager* mcManager = GetOrFatal<CbmMCDataManager>("MCDataManager", description);
CbmMCDataArray* array = mcManager->InitBranch(objName.c_str());
if (array == nullptr) { LOG(fatal) << description << " No MCTrack!"; }
return array;
}
void SaveCanvasAsImage(TCanvas* c, const string& dir, const string& option)
{
if (dir == "") return;
......
......@@ -5,6 +5,11 @@
#ifndef CBMUTILS_H_
#define CBMUTILS_H_
#include "../CbmMCDataManager.h"
#include "FairLogger.h"
#include "FairRootManager.h"
#include <sstream> // for string, stringstream
#include <vector> // for vector
......@@ -46,6 +51,28 @@ namespace Cbm
return (x > ZERO) ? 1 : ((x < ZERO) ? -1 : 0);
}
/*
\brief Tries to get object from FairRootManager. If object is not found create Fatal error.
\param[in] objName Name of the object.
\param[in] description Optional description for LOG, usually class name + method name.
*/
template<typename T>
T* GetOrFatal(const std::string& objName, const std::string& description = "")
{
FairRootManager* ioman = FairRootManager::Instance();
if (ioman == nullptr) { LOG(fatal) << description << " No FairRootManager!"; }
T* obj = static_cast<T*>(ioman->GetObject(objName.c_str()));
if (obj == nullptr) { LOG(fatal) << description << " No " << objName << " object!"; }
return obj;
}
/*
\brief Tries to get CbmMCDataArray from CbmMCDataManager. If object is not found create Fatal error.
\param[in] objName Name of the object.
\param[in] description Optional description for LOG, usually class name + method name.
*/
CbmMCDataArray* InitOrFatalMc(const std::string& objName, const std::string& description = "");
void SaveCanvasAsImage(TCanvas* c, const std::string& dir, const std::string& option = "eps;png;gif");
void SaveCanvasAsImageImpl(const std::string& imageType, TCanvas* c, const std::string& dir,
const std::string& option);
......
......@@ -4,7 +4,9 @@ set(INCLUDE_DIRECTORIES
${CBMDATA_DIR}
${CBMDATA_DIR}/rich
${CBMDATA_DIR}/tof
${CBMDATA_DIR}/global
${CBMDATA_DIR}/base
${CBMBASE_DIR}
${CBMBASE_DIR}/draw
......@@ -39,6 +41,7 @@ set(NO_DICT_SRCS
CbmRichGeoManager.cxx
CbmRichDigiMapManager.cxx
CbmRichMCbmDigiMapManager.cxx
utils/CbmRichUtil.cxx
)
IF (SSE_FOUND)
......
......@@ -5,9 +5,8 @@
#ifndef RICH_CbmRichDraw
#define RICH_CbmRichDraw
#include "CbmRichDetectorData.h" // for CbmRichPmtData, CbmRichPixelData
#include "CbmRichDigiMapManager.h" // for CbmRichDigiMapManager
#include "CbmRichGeoManager.h" // for CbmRichGeoManager
#include "CbmRichUtil.h"
#include <RtypesCore.h> // for ROOT data types
#include <TCanvas.h> // for TCanvas
......@@ -29,7 +28,7 @@ public:
TH2D* hUp = (TH2D*) h->Clone();
DrawH2(hUp);
if (usePmtBins) {
std::vector<Double_t> yPmtBins = CbmRichDraw::GetPmtHistYbins();
std::vector<Double_t> yPmtBins = CbmRichUtil::GetPmtHistYbins();
hUp->GetYaxis()->SetRange(yPmtBins.size() / 2 + 1, yPmtBins.size());
}
else {
......@@ -42,7 +41,7 @@ public:
c->cd(2);
TH2D* hDown = (TH2D*) h->Clone();
if (usePmtBins) {
std::vector<Double_t> yPmtBins = CbmRichDraw::GetPmtHistYbins();
std::vector<Double_t> yPmtBins = CbmRichUtil::GetPmtHistYbins();
hDown->GetYaxis()->SetRange(0, yPmtBins.size() / 2 - 1);
}
else {
......@@ -72,48 +71,6 @@ public:
gPad->SetLeftMargin(0.1);
gPad->SetRightMargin(0.15);
}
static std::vector<Double_t> GetPmtHistXbins() { return CbmRichDraw::GetPmtHistBins(true); }
static std::vector<Double_t> GetPmtHistYbins() { return CbmRichDraw::GetPmtHistBins(false); }
private:
static std::vector<Double_t> GetPmtHistBins(Bool_t isX)
{
std::vector<Double_t> initVec;
std::vector<Int_t> pmts = CbmRichDigiMapManager::GetInstance().GetPmtIds();
for (Int_t pmtId : pmts) {
CbmRichPmtData* pmtData = CbmRichDigiMapManager::GetInstance().GetPmtDataById(pmtId);
TVector3 inPos(pmtData->fX, pmtData->fY, pmtData->fZ);
TVector3 outPos;
CbmRichGeoManager::GetInstance().RotatePoint(&inPos, &outPos);
if (isX) {
initVec.push_back(outPos.X() - 0.5 * pmtData->fWidth);
initVec.push_back(outPos.X() + 0.5 * pmtData->fWidth);
}
else {
initVec.push_back(outPos.Y() - 0.5 * pmtData->fHeight);
initVec.push_back(outPos.Y() + 0.5 * pmtData->fHeight);
}
}
sort(initVec.begin(), initVec.end());
std::vector<Double_t> uniVec;
for (unsigned int i = 0; i < initVec.size(); i++) {
if (i == 0) uniVec.push_back(initVec[i]);
if (initVec[i] - uniVec[uniVec.size() - 1] > 0.000001) uniVec.push_back(initVec[i]);
}
// cout << "uniVec.size():" << uniVec.size() << endl;
// for (int i = 0; i < uniVec.size(); i++) {
// cout << std::setprecision(9);
// cout << fixed;
// cout << uniVec[i] << " " ;
// }
// cout << endl;
return uniVec;
}
};
#endif
/* Copyright (C) 2021 Justus-Liebig-Universitaet Giessen, Giessen
SPDX-License-Identifier: GPL-3.0-only
Authors: Semen Lebedev [committer]*/
#include "CbmRichUtil.h"
#include "CbmDigiManager.h"
#include "CbmGlobalTrack.h"
#include "CbmMCDataArray.h"
#include "CbmMatchRecoToMC.h"
#include "CbmRichDetectorData.h"
#include "CbmRichDigi.h"
#include "CbmRichDigiMapManager.h"
#include "CbmRichGeoManager.h"
#include "CbmRichHit.h"
#include "CbmRichRing.h"
#include "TClonesArray.h"
using namespace std;
map<pair<int, int>, int> CbmRichUtil::CreateNofHitsInRingMap(TClonesArray* richHits, CbmMCDataArray* richPoints,
CbmMCDataArray* mcTracks, CbmDigiManager* digiMan)
{
map<pair<int, int>, int> nofHitsInRingMap;
int nofRichHits = richHits->GetEntriesFast();
for (int iHit = 0; iHit < nofRichHits; iHit++) {
const CbmRichHit* hit = static_cast<CbmRichHit*>(richHits->At(iHit));
if (hit == nullptr || hit->GetRefId() < 0) continue;
const CbmRichDigi* digi = digiMan->Get<CbmRichDigi>(hit->GetRefId());
if (digi == nullptr) continue;
const CbmMatch* digiMatch = digiMan->GetMatch(ECbmModuleId::kRich, hit->GetRefId());
if (digiMatch == nullptr) continue;
int eventId = digiMatch->GetMatchedLink().GetEntry();
vector<pair<int, int>> motherIds =
CbmMatchRecoToMC::GetMcTrackMotherIdsForRichHit(digiMan, hit, richPoints, mcTracks, eventId);
for (size_t i = 0; i < motherIds.size(); i++) {
nofHitsInRingMap[motherIds[i]]++;
}
}
return nofHitsInRingMap;
}
double CbmRichUtil::GetRingTrackDistance(Int_t globalTrackId) { return GetRingTrackDistanceImpl(globalTrackId)[0]; }
double CbmRichUtil::GetRingTrackDistanceX(Int_t globalTrackId) { return GetRingTrackDistanceImpl(globalTrackId)[1]; }
double CbmRichUtil::GetRingTrackDistanceY(Int_t globalTrackId) { return GetRingTrackDistanceImpl(globalTrackId)[2]; }
vector<Double_t> CbmRichUtil::GetRingTrackDistanceImpl(Int_t globalTrackId)
{
vector<Double_t> errorVec = {999., 999., 999.};
FairRootManager* ioman = FairRootManager::Instance();
if (ioman == NULL) return errorVec;
// Do we really need static here, depends on ioman->GetObject() method
static TClonesArray* globalTracks = NULL;
static TClonesArray* richRings = NULL;
static TClonesArray* richProjections = NULL;
if (globalTracks == NULL || richRings == NULL || richProjections == NULL) {
//cout << "globalTracks == NULL || richRings == NULL || richProjections == NULL" << endl;
globalTracks = (TClonesArray*) ioman->GetObject("GlobalTrack");
richRings = (TClonesArray*) ioman->GetObject("RichRing");
richProjections = (TClonesArray*) ioman->GetObject("RichProjection");
}
else {
//cout << "globalTracks, richRings, richProjections NOT NULL" << endl;
}
if (globalTracks == NULL || richRings == NULL || richProjections == NULL) {
LOG(error) << "CbmRichUtil::GetRingTrackDistance globalTracks, "
"richRings, richProjections NOT INITIALIZED"
<< endl;
return errorVec;
}
const CbmGlobalTrack* globalTrack = static_cast<const CbmGlobalTrack*>(globalTracks->At(globalTrackId));
if (globalTrack == NULL) return errorVec;
Int_t stsId = globalTrack->GetStsTrackIndex();
if (stsId < 0) return errorVec;
const FairTrackParam* pTrack = static_cast<const FairTrackParam*>(richProjections->At(stsId));
if (pTrack == NULL) return errorVec;
if (pTrack->GetX() == 0 && pTrack->GetY() == 0) return errorVec;
Int_t richId = globalTrack->GetRichRingIndex();
if (richId < 0) return errorVec;
const CbmRichRing* richRing = static_cast<const CbmRichRing*>(richRings->At(richId));
if (richRing == NULL) return errorVec;
// Double_t xRing = richRing->GetCenterX();
// Double_t yRing = richRing->GetCenterY();
Double_t dx = richRing->GetCenterX() - pTrack->GetX();
Double_t dy = richRing->GetCenterY() - pTrack->GetY();
Double_t dist = TMath::Sqrt(dx * dx + dy * dy);
vector<Double_t> v = {dist, dx, dy};
return v;
}
vector<double> CbmRichUtil::GetPmtHistXbins() { return GetPmtHistBins(true); }
vector<double> CbmRichUtil::GetPmtHistYbins() { return GetPmtHistBins(false); }
vector<double> CbmRichUtil::GetPmtHistBins(bool isX)
{
vector<double> initVec;
vector<int> pmts = CbmRichDigiMapManager::GetInstance().GetPmtIds();
for (int pmtId : pmts) {
CbmRichPmtData* pmtData = CbmRichDigiMapManager::GetInstance().GetPmtDataById(pmtId);
TVector3 inPos(pmtData->fX, pmtData->fY, pmtData->fZ);
TVector3 outPos;
CbmRichGeoManager::GetInstance().RotatePoint(&inPos, &outPos);
if (isX) {
initVec.push_back(outPos.X() - 0.5 * pmtData->fWidth);
initVec.push_back(outPos.X() + 0.5 * pmtData->fWidth);
}
else {
initVec.push_back(outPos.Y() - 0.5 * pmtData->fHeight);
initVec.push_back(outPos.Y() + 0.5 * pmtData->fHeight);
}
}
sort(initVec.begin(), initVec.end());
vector<double> uniVec;
for (size_t i = 0; i < initVec.size(); i++) {
if (i == 0) uniVec.push_back(initVec[i]);
if (initVec[i] - uniVec[uniVec.size() - 1] > 0.000001) uniVec.push_back(initVec[i]);
}
return uniVec;
}
/* Copyright (C) 2017-2019 GSI/JINR-LIT, Darmstadt/Dubna
/* Copyright (C) 2017-2021 GSI/JINR-LIT, Darmstadt/Dubna
SPDX-License-Identifier: GPL-3.0-only
Authors: Andrey Lebedev [committer], Semen Lebedev */
#ifndef RICH_CbmRichUtil
#define RICH_CbmRichUtil
#include "CbmGlobalTrack.h"
#include <map>
#include <vector>
#include "FairRootManager.h"
#include "FairTrackParam.h"
#include <Logger.h>
#include "TCanvas.h"
#include "TClonesArray.h"
#include "TH2.h"
#include "TObject.h"
#include <string>
using namespace std;
class CbmDigiManager;
class TClonesArray;
class CbmMCDataArray;
class TH2D;
class CbmRichUtil {
public:
static Double_t GetRingTrackDistance(Int_t globalTrackId)
{
vector<Double_t> v = GetRingTrackDistanceImpl(globalTrackId);
return v[0];
}
static double GetRingTrackDistance(int globalTrackId);
static double GetRingTrackDistanceX(int globalTrackId);
static double GetRingTrackDistanceY(int globalTrackId);
static Double_t GetRingTrackDistanceX(Int_t globalTrackId)
{
vector<Double_t> v = GetRingTrackDistanceImpl(globalTrackId);
return v[1];
}
// Create PMT XY histograms
static std::vector<double> GetPmtHistXbins();
static std::vector<double> GetPmtHistYbins();
static std::vector<double> GetPmtHistBins(bool isX);
static Double_t GetRingTrackDistanceY(Int_t globalTrackId)
{
vector<Double_t> v = GetRingTrackDistanceImpl(globalTrackId);
return v[2];
}
static std::map<std::pair<int, int>, int> CreateNofHitsInRingMap(TClonesArray* richHits, CbmMCDataArray* richPoints,
CbmMCDataArray* mcTracks, CbmDigiManager* digiMan);
static uint16_t GetDirichId(Int_t Address) { return ((Address >> 16) & 0xFFFF); }
static uint16_t GetDirichId(int Address) { return ((Address >> 16) & 0xFFFF); }
static uint16_t GetDirichChannel(Int_t Address) { return (Address & 0xFFFF); }
static uint16_t GetDirichChannel(int Address) { return (Address & 0xFFFF); }
private:
/**
* \brief Return a vector with total distance and x, y components. [0] - total distance, [1] - x component, [2] - y component
*/
static vector<Double_t> GetRingTrackDistanceImpl(Int_t globalTrackId)
{
vector<Double_t> errorVec = {999., 999., 999.};
FairRootManager* ioman = FairRootManager::Instance();
if (ioman == NULL) return errorVec;
// Do we really need static here, depends on ioman->GetObject() method
static TClonesArray* globalTracks = NULL;
static TClonesArray* richRings = NULL;
static TClonesArray* richProjections = NULL;
if (globalTracks == NULL || richRings == NULL || richProjections == NULL) {
//cout << "globalTracks == NULL || richRings == NULL || richProjections == NULL" << endl;
globalTracks = (TClonesArray*) ioman->GetObject("GlobalTrack");
richRings = (TClonesArray*) ioman->GetObject("RichRing");
richProjections = (TClonesArray*) ioman->GetObject("RichProjection");
}
else {
//cout << "globalTracks, richRings, richProjections NOT NULL" << endl;
}
if (globalTracks == NULL || richRings == NULL || richProjections == NULL) {
LOG(error) << "CbmRichUtil::GetRingTrackDistance globalTracks, "
"richRings, richProjections NOT INITIALIZED"
<< endl;
return errorVec;
}
const CbmGlobalTrack* globalTrack = static_cast<const CbmGlobalTrack*>(globalTracks->At(globalTrackId));
if (globalTrack == NULL) return errorVec;
Int_t stsId = globalTrack->GetStsTrackIndex();
if (stsId < 0) return errorVec;
const FairTrackParam* pTrack = static_cast<const FairTrackParam*>(richProjections->At(stsId));
if (pTrack == NULL) return errorVec;
if (pTrack->GetX() == 0 && pTrack->GetY() == 0) return errorVec;
Int_t richId = globalTrack->GetRichRingIndex();
if (richId < 0) return errorVec;
const CbmRichRing* richRing = static_cast<const CbmRichRing*>(richRings->At(richId));
if (richRing == NULL) return errorVec;
// Double_t xRing = richRing->GetCenterX();
// Double_t yRing = richRing->GetCenterY();
Double_t dx = richRing->GetCenterX() - pTrack->GetX();
Double_t dy = richRing->GetCenterY() - pTrack->GetY();
Double_t dist = TMath::Sqrt(dx * dx + dy * dy);
vector<Double_t> v = {dist, dx, dy};
return v;
}
static std::vector<double> GetRingTrackDistanceImpl(int globalTrackId);
};
#endif
......@@ -11,4 +11,4 @@ commandXterm3 = ('xterm -hold -e python3 run_one.py {} {}&').format(3, "10gev")
os.system(commandXterm3)
commandXterm4 = ('xterm -hold -e python3 run_one.py {} {}&').format(4, "12gev")
os.system(commandXterm4)
\ No newline at end of file
os.system(commandXterm4)
......@@ -32,9 +32,9 @@ def main():
recoCmd = ('root -l -b -q run_reco.C\(\\"{}\\",\\"{}\\",\\"{}\\",\\"{}\\",\\"{}\\",\\"{}\\",{}\)').format(testType, traFile, parFile, digiFile, recoFile, geoSetup, nEvents)
qaCmd = ('root -l -b -q run_qa.C\(\\"{}\\",\\"{}\\",\\"{}\\",\\"{}\\",\\"{}\\",\\"{}\\",\\"{}\\",\\"{}\\",{}\)').format(testType, traFile, parFile, digiFile, recoFile, qaFile, geoSetup, resultDir, nEvents)
os.system((". /{} -a; {}").format(cbmrootConfigPath, traCmd))
os.system((". /{} -a; {}").format(cbmrootConfigPath, digiCmd))
os.system((". /{} -a; {}").format(cbmrootConfigPath, recoCmd))
# os.system((". /{} -a; {}").format(cbmrootConfigPath, traCmd))
# os.system((". /{} -a; {}").format(cbmrootConfigPath, digiCmd))
# os.system((". /{} -a; {}").format(cbmrootConfigPath, recoCmd))
os.system((". /{} -a; {}").format(cbmrootConfigPath, qaCmd))
# def make_args():
......
......@@ -20,6 +20,7 @@
#include "CbmMCTrack.h"
#include "CbmReport.h"
#include "CbmRichConverter.h"
#include "CbmRichDetectorData.h"
#include "CbmRichDigi.h"
#include "CbmRichDigiMapManager.h"
#include "CbmRichDraw.h"
......
This diff is collapsed.
/* Copyright (C) 2012-2020 UGiessen/JINR-LIT, Giessen/Dubna
/* Copyright (C) 2012-2021 UGiessen/JINR-LIT, Giessen/Dubna
SPDX-License-Identifier: GPL-3.0-only
Authors: Andrey Lebedev, Semen Lebedev [committer] */
/**
* \file CbmRichUrqmdTest.h
*
* \brief RICH geometry testing in Urqmd collisions.
*
* \author Semen Lebedev <s.lebedev@gsi.de>
* \date 2012
**/
#ifndef CBM_RICH_URQMD_TEST
#define CBM_RICH_URQMD_TEST
......@@ -26,16 +17,7 @@ class CbmDigiManager;
#include <map>
#include <vector>
using namespace std;
/**
* \class CbmRichUrqmdTest
*
* \brief RICH geometry testing in Urqmd collisions.
*
* \author Semen Lebedev <s.lebedev@gsi.de>
* \date 2012
**/
class CbmRichUrqmdTest : public FairTask {
public:
......@@ -69,20 +51,35 @@ public:
* \brief Set output directory where you want to write results (figures and json).
* \param[in] dir Path to the output directory.
*/
void SetOutputDir(const string& dir) { fOutputDir = dir; }
void SetOutputDir(const std::string& dir) { fOutputDir = dir; }
private:
CbmHistManager* fHM = nullptr;
std::string fOutputDir = ""; // output dir for results
TClonesArray* fRichHits = nullptr;
TClonesArray* fRichRings = nullptr;
CbmMCDataArray* fRichPoints = nullptr;
CbmMCDataArray* fMcTracks = nullptr;
TClonesArray* fRichRingMatches = nullptr;
TClonesArray* fRichProjections = nullptr;
CbmDigiManager* fDigiMan = nullptr;
CbmMCEventList* fEventList = nullptr;
Int_t fEventNum = 0;
Int_t fMinNofHits = 7; // Min number of hits in ring for detector acceptance calculation.
std::map<std::pair<Int_t, Int_t>, Int_t> fNofHitsInRingMap; // Number of hits in the MC RICH ring
std::vector<std::pair<int, int>> fVertexZStsSlices; // Slices (Zmin - Zmax) inside STS
/**
* \brief Initialize histograms.
*/
void InitHistograms();
/**
* \brief
*/
void FillRichRingNofHits();
/**
* \brief
*/
......@@ -121,25 +118,6 @@ private:
*/
CbmRichUrqmdTest& operator=(const CbmRichUrqmdTest&);
CbmHistManager* fHM;
string fOutputDir; // output dir for results
TClonesArray* fRichHits;
TClonesArray* fRichRings;
CbmMCDataArray* fRichPoints;
CbmMCDataArray* fMcTracks;
TClonesArray* fRichRingMatches;
TClonesArray* fRichProjections;
CbmDigiManager* fDigiMan;
CbmMCEventList* fEventList;
Int_t fEventNum;
Int_t fMinNofHits; // Min number of hits in ring for detector acceptance calculation.
// Number of hits in the MC RICH ring
std::map<pair<Int_t, Int_t>, Int_t> fNofHitsInRingMap;
ClassDef(CbmRichUrqmdTest, 1)
};
......