-
Sergey Gorbunov authoredSergey Gorbunov authored
CbmL1.cxx 43.94 KiB
/* Copyright (C) 2006-2023 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
SPDX-License-Identifier: GPL-3.0-only
Authors: Ivan Kisel, Sergey Gorbunov [committer], Valentina Akishina, Igor Kulakov, Maksym Zyzak, Sergei Zharko */
/*
*====================================================================
*
* CBM Level 1 Reconstruction
*
* Authors: I.Kisel, S.Gorbunov
*
* e-mail : ikisel@kip.uni-heidelberg.de
*
*====================================================================
*
* L1 main program
*
*====================================================================
*/
#include "CbmL1.h"
#include "CaToolsMaterialHelper.h"
#include "CbmMCDataManager.h"
#include "CbmMuchTrackingInterface.h"
#include "CbmMvdTrackingInterface.h"
#include "CbmSetup.h"
#include "CbmStsTrackingInterface.h"
#include "CbmTofTrackingInterface.h"
#include "CbmTrdTrackingInterface.h"
#include "KfDefs.h"
#include <boost/filesystem.hpp>
// TODO: include of CbmSetup.h creates problems on Mac
// #include "CbmSetup.h"
#include "CaFramework.h"
#include "CaHit.h"
#include "CaToolsDebugger.h"
#include "CaToolsField.h"
#include "CbmEvent.h"
#include "CbmKfTrackingSetupInitializer.h"
#include "CbmMCDataObject.h"
#include "CbmStsFindTracks.h"
#include "CbmStsHit.h"
#include "CbmTrackingDetectorInterfaceInit.h"
#include "FairField.h"
#include "FairRunAna.h"
#include "TGeoArb8.h"
#include "TGeoBoolNode.h"
#include "TGeoCompositeShape.h"
#include "TGeoManager.h"
#include "TGeoNode.h"
#include "TMatrixD.h"
#include "TNtuple.h"
#include "TProfile2D.h"
#include "TROOT.h"
#include "TRandom3.h"
#include "TSystem.h"
#include "TVector3.h"
#include "TVectorD.h"
#include <TFile.h>
#include <chrono>
#include <fstream>
#include <iomanip>
#include <list>
#include <sstream>
using namespace cbm::algo; // TODO: remove this line
using CaTrack = cbm::algo::ca::Track;
using cbm::algo::ca::Parameters;
using cbm::ca::MCModule;
using cbm::ca::TimeSliceReader;
using cbm::ca::tools::MaterialHelper;
using cbm::kf::TrackingSetupInitializer;
ClassImp(CbmL1);
static ca::Framework gAlgo _fvecalignment; // TODO: Change coupling logic between ca::Framework and CbmL1
CbmL1* CbmL1::fpInstance = nullptr;
// ---------------------------------------------------------------------------------------------------------------------
//
CbmL1::CbmL1() : CbmL1("L1") {}
// ---------------------------------------------------------------------------------------------------------------------
//
CbmL1::CbmL1(const char* name, Int_t verbose, Int_t performance) : FairTask(name, verbose), fPerformance(performance)
{
if (!fpInstance) fpInstance = this;
switch (fSTAPDataMode) {
case 1: LOG(info) << "CbmL1: input data will be written for a standalone usage"; break;
case 2: LOG(info) << "CbmL1: input data will be read from external files"; break;
default: LOG(info) << "CbmL1: tracking will be run without external data R/W"; break;
}
fpIODataManager = std::make_shared<ca::DataManager>();
this->DefineSTAPNames("");
if (!CbmTrackingDetectorInterfaceInit::Instance()) {
LOG(fatal) << "CbmL1: CbmTrackingDetectorInterfaceInit instance was not found. Please, add it as a task to your "
"reco macro right before the KF and L1 tasks:\n"
<< "\033[1;30mrun->AddTask(new CbmTrackingDetectorInterfaceInit());\033[0m";
}
}
// ---------------------------------------------------------------------------------------------------------------------
//
CbmL1::~CbmL1()
{
if (fpInstance == this) fpInstance = nullptr;
}
// ---------------------------------------------------------------------------------------------------------------------
//
void CbmL1::CheckDetectorPresence()
{
fUseMVD = fUseMVD && CbmSetup::Instance()->IsActive(ECbmModuleId::kMvd);
fUseSTS = fUseSTS && CbmSetup::Instance()->IsActive(ECbmModuleId::kSts);
fUseMUCH = fUseMUCH && CbmSetup::Instance()->IsActive(ECbmModuleId::kMuch);
fUseTRD = fUseTRD && CbmSetup::Instance()->IsActive(ECbmModuleId::kTrd);
fUseTOF = fUseTOF && CbmSetup::Instance()->IsActive(ECbmModuleId::kTof);
{
// TODO: temporary code!!
// for a moment, the MVD digitizer doesn't work in TB mode
// check the presence of MVD hits to make sure the MVD is really active
if (!FairRootManager::Instance()->GetObject("MvdHit")) {
fUseMVD = false;
}
}
}
// ---------------------------------------------------------------------------------------------------------------------
//void CbmL1::DisableTrackingStation(ca::EDetectorID detID, int iSt)
{
if (ca::EDetectorID::kEND != detID) {
fvmDisabledStationIDs[detID].insert(iSt);
}
}
// ---------------------------------------------------------------------------------------------------------------------
//
InitStatus CbmL1::ReInit()
{
SetParContainers();
return Init();
}
// ---------------------------------------------------------------------------------------------------------------------
//
InitStatus CbmL1::Init()
try {
if (fVerbose > 1) {
char y[20] = " [0;33;44m"; // yellow
char Y[20] = " [1;33;44m"; // yellow bold
char W[20] = " [1;37;44m"; // white bold
char o[20] = " [0m\n"; // color off
Y[0] = y[0] = W[0] = o[0] = 0x1B; // escape character
std::stringstream ss;
ss << "\n\n";
ss << " " << W << " " << o;
ss << " " << W << " ===////====================================================== " << o;
ss << " " << W << " = = " << o;
ss << " " << W << " = " << Y << "L1 on-line finder" << W << " = " << o;
ss << " " << W << " = = " << o;
ss << " " << W << " = " << W << "Cellular Automaton 3.1 Vector" << y << " with " << W << "KF Quadro" << y
<< " technology" << W << " = " << o;
ss << " " << W << " = = " << o;
ss << " " << W << " = " << y << "Designed for CBM collaboration" << W << " = " << o;
ss << " " << W << " = " << y << "All rights reserved" << W << " = " << o;
ss << " " << W << " = = " << o;
ss << " " << W << " ========================================================////= " << o;
ss << " " << W << " " << o;
ss << "\n\n";
LOG(info) << ss.str();
}
cbm::ca::tools::SetOriginalCbmField();
fHistoDir = gROOT->mkdir("L1");
fHistoDir->mkdir("Input");
fHistoDir->mkdir("Fit");
fTableDir = gROOT->mkdir("L1Tables");
// turn on reconstruction in sub-detectors
fUseMVD = false;
fUseSTS = false;
fUseMUCH = false;
fUseTRD = false;
fUseTOF = false;
if (ca::TrackingMode::kSts == fTrackingMode) {
fUseMVD = true;
fUseSTS = true;
fUseMUCH = false;
fUseTRD = false;
fUseTOF = false;
// check if MVD is switched off in the Sts task
CbmStsFindTracks* findTask = dynamic_cast<CbmStsFindTracks*>(FairRunAna::Instance()->GetTask("STSFindTracks")); if (findTask) fUseMVD = findTask->MvdUsage();
}
if (ca::TrackingMode::kMcbm == fTrackingMode) {
fUseMVD = false;
fUseSTS = true;
fUseMUCH = true;
fUseTRD = true;
fUseTOF = true;
// fInitManager.DevSetIgnoreHitSearchAreas(true); // uncomment for debug
}
if (ca::TrackingMode::kGlobal == fTrackingMode) {
//at the moment trd2d tracking only
fUseMVD = false;
fUseSTS = false;
fUseMUCH = false;
fUseTRD = true;
fUseTOF = false;
fInitManager.DevSetUseOfOriginalField();
fInitManager.DevSetIgnoreHitSearchAreas(true);
fInitManager.SetMaxTripletPerDoublets(1000);
}
// uncomment for debug
//fInitManager.DevSetIgnoreHitSearchAreas(true);
//fInitManager.DevSetIsMatchDoubletsViaMc(true);
//fInitManager.DevSetIsMatchTripletsViaMc(true);
//fInitManager.DevSetIsExtendTracksViaMc(true);
//fInitManager.DevSetIsSuppressOverlapHitsViaMc(true);
CheckDetectorPresence();
// *****************************
// ** **
// ** GEOMETRY INITIALIZATION **
// ** **
// *****************************
// Read parameters object from a binary file
if (2 == fSTAPDataMode) {
this->ReadSTAPParamObject();
}
// Parameters initialization, based on the FairRunAna chain
else {
// ********************************************
// ** Base configuration file initialization **
// ********************************************
{
// The parameters module of CA algorithm (L1Parameters) contains two sets of data fields: the geometry configura-
// tion and the parameter configuration. The configuration both for geometry and parameters can be set using
// either the accessor functions of the L1InitManager class, or reading the beforehand prepared L1Parameters
// object. On top of that, the parameters configuration (including cuts, algorithm execution scenario etc.) can
// be initialized from the YAML configuration files in two levels (base and user). The base level of the initia-
// lization is required for CbmL1 mode, the user level of the initialization is optional.
std::string mainConfig = std::string(gSystem->Getenv("VMCWORKDIR")) + "/macro/L1/configs/";
switch (fTrackingMode) {
case ca::TrackingMode::kSts: mainConfig += "ca_params_sts.yaml"; break;
case ca::TrackingMode::kMcbm: mainConfig += "ca_params_mcbm.yaml"; break;
case ca::TrackingMode::kGlobal: mainConfig += "ca_params_global.yaml"; break;
}
fInitManager.SetConfigMain(mainConfig);
}
fInitManager.SetDetectorNames(cbm::ca::kDetName);
auto mvdInterface = CbmMvdTrackingInterface::Instance(); auto stsInterface = CbmStsTrackingInterface::Instance();
auto muchInterface = CbmMuchTrackingInterface::Instance();
auto trdInterface = CbmTrdTrackingInterface::Instance();
auto tofInterface = CbmTofTrackingInterface::Instance();
int nMvdStationsGeom = (fUseMVD) ? mvdInterface->GetNtrackingStations() : 0;
int nStsStationsGeom = (fUseSTS) ? stsInterface->GetNtrackingStations() : 0;
int nMuchStationsGeom = (fUseMUCH) ? muchInterface->GetNtrackingStations() : 0;
int nTrdStationsGeom = (fUseTRD) ? trdInterface->GetNtrackingStations() : 0;
int nTofStationsGeom = (fUseTOF) ? tofInterface->GetNtrackingStations() : 0;
// **************************
// ** Field initialization **
// **************************
if (FairRunAna::Instance()->GetField()) {
fInitManager.SetFieldFunction([](const double(&inPos)[3], double(&outB)[3]) {
FairRunAna::Instance()->GetField()->GetFieldValue(inPos, outB);
});
}
else {
fInitManager.SetFieldFunction([](const double(&)[3], double(&outB)[3]) {
outB[0] = 0.;
outB[1] = 0.;
outB[2] = 0.;
});
}
// ***************************
// ** Target initialization **
// ***************************
GetTargetInfo();
fInitManager.SetTargetPosition(fTargetX, fTargetY, fTargetZ);
// *********************************
// ** Target field initialization **
// *********************************
// FIXME: SZh 22.08.2023:
// Is there anyway to calculate a step between field slices?
fInitManager.InitTargetField(/*zStep = */ 2.5 /*cm*/); // Replace zStep -> sizeZfieldRegion = 2 * zStep (TODO)
// **************************************
// ** **
// ** STATIONS GEOMETRY INITIALIZATION **
// ** **
// **************************************
// ***************************************************
// ** Active tracking detector subsystems selection **
// ***************************************************
std::set<ca::EDetectorID> vActiveTrackingDetectorIDs{}; // Set of detectors active in tracking
if (fUseMVD) {
vActiveTrackingDetectorIDs.insert(ca::EDetectorID::kMvd);
}
if (fUseSTS) {
vActiveTrackingDetectorIDs.insert(ca::EDetectorID::kSts);
}
if (fUseMUCH) {
vActiveTrackingDetectorIDs.insert(ca::EDetectorID::kMuch);
}
if (fUseTRD) {
vActiveTrackingDetectorIDs.insert(ca::EDetectorID::kTrd);
} if (fUseTOF) {
vActiveTrackingDetectorIDs.insert(ca::EDetectorID::kTof);
}
fInitManager.SetActiveDetectorIDs(vActiveTrackingDetectorIDs);
// *************************************
// ** Stations layout initialization **
// *************************************
std::vector<ca::StationInitializer> vStations;
vStations.reserve(100);
bool bDisableTime = false; /// TMP, move to parameters!
// *** MVD stations info ***
if (fUseMVD) {
for (int iSt = 0; iSt < nMvdStationsGeom; ++iSt) {
auto stationInfo = ca::StationInitializer(ca::EDetectorID::kMvd, iSt);
// TODO: SZh 15.08.2022: Replace station type with ca::EDetectorID
stationInfo.SetStationType(1); // MVD
stationInfo.SetTimeInfo(!bDisableTime && mvdInterface->IsTimeInfoProvided(iSt));
stationInfo.SetFieldStatus(fTrackingMode == ca::TrackingMode::kMcbm ? 0 : 1);
stationInfo.SetZref(mvdInterface->GetZref(iSt));
stationInfo.SetZmin(mvdInterface->GetZmin(iSt));
stationInfo.SetZmax(mvdInterface->GetZmax(iSt));
stationInfo.SetXmax(std::max(fabs(mvdInterface->GetXmax(iSt)), fabs(mvdInterface->GetXmin(iSt))));
stationInfo.SetYmax(std::max(fabs(mvdInterface->GetYmax(iSt)), fabs(mvdInterface->GetYmin(iSt))));
stationInfo.SetTrackingStatus(true);
if (fvmDisabledStationIDs[ca::EDetectorID::kMvd].find(iSt)
!= fvmDisabledStationIDs[ca::EDetectorID::kMvd].cend()) {
stationInfo.SetTrackingStatus(false);
}
fInitManager.AddStation(stationInfo);
LOG(info) << "- MVD station " << iSt << " at z = " << stationInfo.GetZref() << " cm ";
}
}
// *** STS stations info ***
if (fUseSTS) {
for (int iSt = 0; iSt < nStsStationsGeom; ++iSt) {
auto stationInfo = ca::StationInitializer(ca::EDetectorID::kSts, iSt);
// TODO: SZh 15.08.2022: Replace station type with ca::EDetectorID
stationInfo.SetStationType(0); // STS
stationInfo.SetTimeInfo(!bDisableTime && stsInterface->IsTimeInfoProvided(iSt));
stationInfo.SetFieldStatus(ca::TrackingMode::kMcbm == fTrackingMode ? 0 : 1);
stationInfo.SetZref(stsInterface->GetZref(iSt));
stationInfo.SetZmin(stsInterface->GetZmin(iSt));
stationInfo.SetZmax(stsInterface->GetZmax(iSt));
stationInfo.SetXmax(std::max(fabs(stsInterface->GetXmax(iSt)), fabs(stsInterface->GetXmin(iSt))));
stationInfo.SetYmax(std::max(fabs(stsInterface->GetYmax(iSt)), fabs(stsInterface->GetYmin(iSt))));
stationInfo.SetTrackingStatus(true);
if (fvmDisabledStationIDs[ca::EDetectorID::kSts].find(iSt)
!= fvmDisabledStationIDs[ca::EDetectorID::kSts].cend()) {
stationInfo.SetTrackingStatus(false);
}
fInitManager.AddStation(stationInfo);
LOG(info) << "- STS station " << iSt << " at z = " << stationInfo.GetZref() << " cm ";
}
}
// *** MuCh stations info ***
if (fUseMUCH) {
for (int iSt = 0; iSt < nMuchStationsGeom; ++iSt) {
auto stationInfo = ca::StationInitializer(ca::EDetectorID::kMuch, iSt);
// TODO: SZh 15.08.2022: Replace station type with ca::EDetectorID
stationInfo.SetStationType(2); // MuCh
stationInfo.SetTimeInfo(!bDisableTime && muchInterface->IsTimeInfoProvided(iSt));
stationInfo.SetFieldStatus(0);
stationInfo.SetZref(muchInterface->GetZref(iSt)); stationInfo.SetZmin(muchInterface->GetZmin(iSt));
stationInfo.SetZmax(muchInterface->GetZmax(iSt));
stationInfo.SetXmax(std::max(fabs(muchInterface->GetXmax(iSt)), fabs(muchInterface->GetXmin(iSt))));
stationInfo.SetYmax(std::max(fabs(muchInterface->GetYmax(iSt)), fabs(muchInterface->GetYmin(iSt))));
stationInfo.SetTrackingStatus(true);
if (fvmDisabledStationIDs[ca::EDetectorID::kMuch].find(iSt)
!= fvmDisabledStationIDs[ca::EDetectorID::kMuch].cend()) {
stationInfo.SetTrackingStatus(false);
}
fInitManager.AddStation(stationInfo);
LOG(info) << "- MuCh station " << iSt << " at z = " << stationInfo.GetZref() << " cm";
}
}
// *** TRD stations info ***
if (fUseTRD) {
for (int iSt = 0; iSt < nTrdStationsGeom; ++iSt) {
auto stationInfo = ca::StationInitializer(ca::EDetectorID::kTrd, iSt);
// TODO: SZh 15.08.2022: Replace station type with ca::EDetectorID
stationInfo.SetStationType(3);
stationInfo.SetTimeInfo(!bDisableTime && trdInterface->IsTimeInfoProvided(iSt));
stationInfo.SetFieldStatus(0);
stationInfo.SetZref(trdInterface->GetZref(iSt));
stationInfo.SetZmin(trdInterface->GetZmin(iSt));
stationInfo.SetZmax(trdInterface->GetZmax(iSt));
stationInfo.SetXmax(std::max(fabs(trdInterface->GetXmax(iSt)), fabs(trdInterface->GetXmin(iSt))));
stationInfo.SetYmax(std::max(fabs(trdInterface->GetYmax(iSt)), fabs(trdInterface->GetYmin(iSt))));
if (ca::TrackingMode::kGlobal == fTrackingMode) {
stationInfo.SetTimeInfo(false);
}
stationInfo.SetTrackingStatus(true);
if (fvmDisabledStationIDs[ca::EDetectorID::kTrd].find(iSt)
!= fvmDisabledStationIDs[ca::EDetectorID::kTrd].cend()) {
stationInfo.SetTrackingStatus(false);
}
fInitManager.AddStation(stationInfo);
LOG(info) << "- TRD station " << iSt << " at z = " << stationInfo.GetZref() << " cm";
}
}
// *** TOF stations info ***
if (fUseTOF) {
for (int iSt = 0; iSt < nTofStationsGeom; ++iSt) {
auto stationInfo = ca::StationInitializer(ca::EDetectorID::kTof, iSt);
// TODO: SZh 15.08.2022: Replace station type with ca::EDetectorID
stationInfo.SetStationType(4);
stationInfo.SetTimeInfo(!bDisableTime && tofInterface->IsTimeInfoProvided(iSt));
stationInfo.SetFieldStatus(0);
stationInfo.SetZref(tofInterface->GetZref(iSt));
stationInfo.SetZmin(tofInterface->GetZmin(iSt));
stationInfo.SetZmax(tofInterface->GetZmax(iSt));
stationInfo.SetXmax(std::max(fabs(tofInterface->GetXmax(iSt)), fabs(tofInterface->GetXmin(iSt))));
stationInfo.SetYmax(std::max(fabs(tofInterface->GetYmax(iSt)), fabs(tofInterface->GetYmin(iSt))));
stationInfo.SetTrackingStatus(true);
if (fvmDisabledStationIDs[ca::EDetectorID::kTof].find(iSt)
!= fvmDisabledStationIDs[ca::EDetectorID::kTof].cend()) {
stationInfo.SetTrackingStatus(false);
}
fInitManager.AddStation(stationInfo);
LOG(info) << "- TOF station " << iSt << " at z = " << stationInfo.GetZref() << " cm";
}
}
// Init station layout: sort stations in z-position and init maps of station local, geo and active indices
fInitManager.InitStationLayout();
fInitManager.ReadInputConfigs();
// ****************************************
// ** Material maps initialization **
// ****************************************
this->GenerateMaterialMaps();
// *******************************
// ** Initialize search windows **
// *******************************
if (fsInputSearchWindowsFilename.size()) {
fInitManager.DevSetIsParSearchWUsed();
fInitManager.ReadSearchWindows(fsInputSearchWindowsFilename);
}
else {
fInitManager.DevSetIsParSearchWUsed(false);
}
// Form parameters container
if (!fInitManager.FormParametersContainer()) {
return kFATAL;
}
// Write parameters object to file if needed
if (1 == fSTAPDataMode || 4 == fSTAPDataMode) {
this->WriteSTAPParamObject();
}
}
if (fInitMode == EInitMode::Param) {
auto parameters = fInitManager.TakeParameters();
LOG(info) << '\n' << parameters.ToString(1);
return kSUCCESS;
}
// Init L1 algo core
// FIXME: SZh 22.08.2023:
// Re-organize the the relation between CbmL1 and ca::Framework. I believe, we don't need a global pointer to the ca::Framework
// instance.
fpAlgo = &gAlgo;
//
// ** Send formed parameters object to ca::Framework instance **
//
{
auto parameters = fInitManager.TakeParameters();
fpParameters = std::make_shared<ca::Parameters<ca::fvec>>(parameters);
fpAlgo->ReceiveParameters(std::move(parameters));
}
fpAlgo->Init(fTrackingMode);
if constexpr (0) {
// **** FEATURE: KF-SETUP INITIALIZATION *****
// Creating a setup
auto pSetupInitializer = std::make_unique<TrackingSetupInitializer>();
pSetupInitializer->Use(fUseMVD, fUseSTS, fUseMUCH, fUseTRD, fUseTOF);
pSetupInitializer->InitProperties();
auto trackerSetup = pSetupInitializer->MakeSetup<double>();
LOG(info) << "!!!!\n!!!!\n!!!!\n!!!! KF SETUP:\n" << trackerSetup.ToString(0) << "\n!!!!\n!!!!\n!!!!";
// Storing setup
TrackingSetupInitializer::Store(trackerSetup, "./trackerSetup.geo.kf.bin");
// Reading setup (now in fvec):
auto trackerSetupVec = TrackingSetupInitializer::Load<fvec>("./trackerSetup.geo.kf.bin");
LOG(info) << "!!!!\n!!!!\n!!!!\n!!!! KF SETUP (from file):\n"
<< trackerSetupVec.ToString(1) << "\n!!!!\n!!!!\n!!!!";
{ // TEST: magnetic field
const auto& fld = trackerSetupVec.GetField();
struct Node {
int iSt;
fvec x;
fvec y;
};
{
Node n{.iSt = 1, .x = 2.3, .y = 1.2};
auto kfFldValue = fld.GetFieldValue(n.iSt, n.x, n.y);
auto caFldValue = fpParameters->GetStation(n.iSt).fieldSlice.GetFieldValue(n.x, n.y);
LOG(info) << "A: " << kfFldValue.ToString() << " --- " << caFldValue.ToString();
}
{
Node n{.iSt = 2, .x = -2.1, .y = 2.2};
auto kfFldValue = fld.GetFieldValue(n.iSt, n.x, n.y);
auto caFldValue = fpParameters->GetStation(n.iSt).fieldSlice.GetFieldValue(n.x, n.y);
LOG(info) << "B: " << kfFldValue.ToString() << " --- " << caFldValue.ToString();
}
{
Node n{.iSt = 5, .x = -2.1, .y = 2.2};
auto kfFldValue = fld.GetFieldValue(n.iSt, n.x, n.y);
auto caFldValue = fpParameters->GetStation(n.iSt).fieldSlice.GetFieldValue(n.x, n.y);
LOG(info) << "B: " << kfFldValue.ToString() << " --- " << caFldValue.ToString();
}
// TEST: field region building
{
using cbm::algo::kf::EFieldType;
fvec x = 1.2;
fvec y = 1.3;
fvec z0 = -32.;
fvec z1 = -28.;
fvec z2 = -24.;
auto kfFldVal_0 = fld.GetFieldValue(0, x, y);
auto kfFldVal_1 = fld.GetFieldValue(1, x, y);
auto kfFldVal_2 = fld.GetFieldValue(2, x, y);
auto kfFldReg = fld.GetFieldRegion(kfFldVal_0, z0, kfFldVal_1, z1, kfFldVal_2, z2);
auto caFldVal_0 = fpParameters->GetStation(0).fieldSlice.GetFieldValue(x, y);
auto caFldVal_1 = fpParameters->GetStation(1).fieldSlice.GetFieldValue(x, y);
auto caFldVal_2 = fpParameters->GetStation(2).fieldSlice.GetFieldValue(x, y);
auto caFldReg = cbm::algo::kf::FieldRegion<fvec>{};
caFldReg.Set(caFldVal_0, z0, caFldVal_1, z1, caFldVal_2, z2);
LOG(info) << "CA ----> " << caFldReg.ToString();
LOG(info) << "KF ----> " << kfFldReg.ToString();
}
// TEST: material maps
{
Node n{.iSt = 1, .x = 2.3, .y = 1.2};
LOG(info) << trackerSetupVec.GetMaterial(n.iSt).GetThickness(n.x, n.y) << " -- "
<< fpParameters->GetMaterialThickness(n.iSt, n.x, n.y);
}
{
Node n{.iSt = 1, .x = -2.3, .y = -1.2};
LOG(info) << trackerSetupVec.GetMaterial(n.iSt).GetThickness(n.x, n.y) << " -- "
<< fpParameters->GetMaterialThickness(n.iSt, n.x, n.y);
}
{
Node n{.iSt = 1, .x = 2.5, .y = 0.2};
LOG(info) << trackerSetupVec.GetMaterial(n.iSt).GetThickness(n.x, n.y) << " -- "
<< fpParameters->GetMaterialThickness(n.iSt, n.x, n.y);
}
{
Node n{.iSt = 4, .x = 5.3, .y = -2.2};
LOG(info) << trackerSetupVec.GetMaterial(n.iSt).GetThickness(n.x, n.y) << " -- "
<< fpParameters->GetMaterialThickness(n.iSt, n.x, n.y);
}
{ Node n{.iSt = 4, .x = 0, .y = -1.2};
LOG(info) << trackerSetupVec.GetMaterial(n.iSt).GetThickness(n.x, n.y) << " -- "
<< fpParameters->GetMaterialThickness(n.iSt, n.x, n.y);
}
}
}
// Initialize time-slice reader
fpTSReader = std::make_unique<TimeSliceReader>();
fpTSReader->SetDetector(ca::EDetectorID::kMvd, fUseMVD);
fpTSReader->SetDetector(ca::EDetectorID::kSts, fUseSTS);
fpTSReader->SetDetector(ca::EDetectorID::kMuch, fUseMUCH);
fpTSReader->SetDetector(ca::EDetectorID::kTrd, fUseTRD);
fpTSReader->SetDetector(ca::EDetectorID::kTof, fUseTOF);
fpTSReader->RegisterParameters(fpParameters);
fpTSReader->RegisterHitIndexContainer(fvExternalHits);
fpTSReader->RegisterIODataManager(fpIODataManager);
fpTSReader->RegisterQaHitContainer(fvHitDebugInfo);
if (!fpTSReader->InitRun()) {
return kFATAL;
}
if (fPerformance) {
fpMCModule = std::make_unique<MCModule>(fVerbose, fPerformance);
fpMCModule->SetDetector(ca::EDetectorID::kMvd, fUseMVD);
fpMCModule->SetDetector(ca::EDetectorID::kSts, fUseSTS);
fpMCModule->SetDetector(ca::EDetectorID::kMuch, fUseMUCH);
fpMCModule->SetDetector(ca::EDetectorID::kTrd, fUseTRD);
fpMCModule->SetDetector(ca::EDetectorID::kTof, fUseTOF);
fpMCModule->RegisterMCData(fMCData);
fpMCModule->RegisterRecoTrackContainer(fvRecoTracks);
fpMCModule->RegisterHitIndexContainer(fvExternalHits);
fpMCModule->RegisterParameters(fpParameters);
fpMCModule->RegisterQaHitContainer(fvHitDebugInfo);
fpMCModule->RegisterFirstHitIndexes(fpTSReader->GetHitFirstIndexDet());
if (!fpMCModule->InitRun()) {
return kFATAL;
}
}
LOG(info) << fpParameters->ToString(1);
LOG(info) << "----- Numbers of stations active in tracking -----";
LOG(info) << " MVD: " << fpParameters->GetNstationsActive(ca::EDetectorID::kMvd);
LOG(info) << " STS: " << fpParameters->GetNstationsActive(ca::EDetectorID::kSts);
LOG(info) << " MuCh: " << fpParameters->GetNstationsActive(ca::EDetectorID::kMuch);
LOG(info) << " TRD: " << fpParameters->GetNstationsActive(ca::EDetectorID::kTrd);
LOG(info) << " TOF: " << fpParameters->GetNstationsActive(ca::EDetectorID::kTof);
LOG(info) << " Total: " << fpParameters->GetNstationsActive();
fNStations = fpParameters->GetNstationsActive();
// monitor the material
{
LOG(info) << "\033[31;1m-------------------- L1 material -----------------------------\033[0m";
fMaterialMonitor.clear();
for (int i = 0; i < fNStations; i++) {
ca::MaterialMonitor m(&(fpAlgo->GetParameters().GetThicknessMaps()[i]), Form("station %d", i));
LOG(info) << m.ToString();
fMaterialMonitor.push_back(m);
}
LOG(info) << "\033[31;1m-------------------------------------------------------------\033[0m";
}
DumpMaterialToFile("L1material.root");
// Initialize monitor
fMonitor.Reset();
return kSUCCESS;
}
catch (const std::exception& err) {
LOG(error) << "CbmL1: initialization failed. Reason: " << err.what();
return kFATAL;
}
void CbmL1::Reconstruct(CbmEvent* event)
{
ca::TrackingMonitorData monitorData;
fpTSReader->ReadEvent(event);
if (fPerformance) {
fpMCModule->InitEvent(event); // reads points and MC tracks
fpMCModule->MatchHits(); // matches points with hits
}
fpAlgo->ReceiveInputData(fpIODataManager->TakeInputData());
// Material monitoring: mark active areas
{
for (const auto& hit : fpAlgo->GetInputData().GetHits()) {
fMaterialMonitor[hit.Station()].MarkActiveBin(hit.X(), hit.Y());
}
}
//LOG(info) << "CHECK: hit ids = " << fvExternalHits.size() << ", hits = " << fpIODataManager->GetNofHits()
//<< ", dbg hits = " << fvHitDebugInfo.size();
fpAlgo->SetMonitorData(monitorData);
if (nullptr != event) {
LOG_IF(debug, fVerbose > 0) << "\n======= Ca Track finder: process event " << event->GetNumber() << " ...";
}
else {
LOG_IF(debug, fVerbose > 0) << "\n======= Ca Track finder: process timeslice ...";
}
fpAlgo->FindTracks();
// IdealTrackFinder();
fTrackingTime = fpAlgo->fCaRecoTime; // TODO: remove (not used)
LOG_IF(debug, fVerbose > 0) << "Ca Track Finder finished, found " << fpAlgo->fRecoTracks.size() << " tracks";
// Update monitor data after the actual tracking
monitorData = fpAlgo->GetMonitorData();
// save reconstructed tracks
fvRecoTracks.clear();
fvRecoTracks.reserve(fpAlgo->fRecoTracks.size());
int trackFirstHit = 0;
// FIXME: Rewrite
for (const auto& caTrk : fpAlgo->fRecoTracks) {
CbmL1Track t;
t.Set(caTrk.fParFirst);
t.TLast.Set(caTrk.fParLast);
t.Tpv.Set(caTrk.fParPV);
t.Hits.clear();
for (int i = 0; i < caTrk.fNofHits; i++) {
int caHitId = fpAlgo->fRecoHits[trackFirstHit + i];
int cbmHitID = fpAlgo->GetInputData().GetHit(caHitId).Id();
t.Hits.push_back(cbmHitID);
}
fvRecoTracks.push_back(t);
trackFirstHit += caTrk.fNofHits;
//fMonitor.IncrementCounter(EMonitorKey::kRecoHit, it->fNofHits); }
//fMonitor.IncrementCounter(EMonitorKey::kRecoTrack, fvRecoTracks.size());
LOG(debug) << "CA Track Finder: " << fpAlgo->fCaRecoTime << " s/sub-ts";
// output performance
if (fPerformance) {
LOG_IF(info, fVerbose) << "Performance...";
fpMCModule->MatchTracks(); // matches reco and MC tracks, fills the MC-truth fields of reco tracks
//
// tracker input performance is moved to external QA tasks.
// InputPerformance() method is not optimised to run with the event builder
// TODO: verify QA tasks and remove InputPerformance()
// InputPerformance();
//
EfficienciesPerformance();
HistoPerformance();
TrackFitPerformance();
if (fsMcTripletsOutputFilename.size()) {
DumpMCTripletsToTree();
}
// TimeHist();
/// WriteSIMDKFData();
LOG_IF(info, fVerbose > 1) << "Tracking performance... done";
}
LOG_IF(debug, fVerbose > 1) << "End of CA";
++fEventNo;
fMonitor.AddMonitorData(monitorData);
}
// ----- Finish CbmStsFitPerformanceTask task -----------------------------
void CbmL1::Finish()
{
if (fPerformance) {
// FieldApproxCheck();
// FieldIntegralCheck();
EfficienciesPerformance(kTRUE);
}
// monitor the material
LOG(info) << "\033[31;1m ***************************\033[0m";
LOG(info) << "\033[31;1m ** CA Tracking monitor **\033[0m";
LOG(info) << "\033[31;1m ***************************\033[0m";
LOG(info) << "\033[31;1m ----- Material budget -------------------------------------- \033[0m";
for (int i = 0; i < fNStations; i++) {
LOG(info) << fMaterialMonitor[i].ToString();
}
LOG(info) << "\033[31;1m -------------------------------------------------------------\033[0m";
// monitor of the reconstructed tracks
LOG(info) << fMonitor.ToString();
TDirectory* curr = gDirectory;
TFile* currentFile = gFile;
// Open output file and write histograms
boost::filesystem::path p = (FairRunAna::Instance()->GetUserOutputFileName()).Data();
std::string dir = p.parent_path().string();
if (dir.empty()) dir = ".";
{
std::string histoOutName = dir + "/L1_histo_" + p.filename().string();
LOG(info) << "\033[31;1mL1 performance histograms will be saved to: \033[0m" << histoOutName; TFile* outfile = new TFile(histoOutName.c_str(), "RECREATE");
outfile->cd();
writedir2current(fHistoDir);
outfile->Close();
outfile->Delete();
}
{
std::string tablesOutName = dir + "/L1_perftable_" + p.filename().string();
LOG(info) << "\033[31;1mL1 performance tables will be saved to: \033[0m" << tablesOutName;
TFile* outfile = new TFile(tablesOutName.c_str(), "RECREATE");
outfile->cd();
writedir2current(fTableDir);
outfile->Close();
outfile->Delete();
}
if (fpMcTripletsOutFile) {
fpMcTripletsOutFile->cd();
fpMcTripletsTree->Write();
fpMcTripletsOutFile->Close();
fpMcTripletsOutFile->Delete();
}
gFile = currentFile;
gDirectory = curr;
fpAlgo->Finish();
cbm::ca::tools::Debugger::Instance().Write();
}
// ---------------------------------------------------------------------------------------------------------------------
//
void CbmL1::writedir2current(TObject* obj)
{
if (!obj->IsFolder())
obj->Write();
else {
TDirectory* cur = gDirectory;
TDirectory* sub = cur->mkdir(obj->GetName());
sub->cd();
TList* listSub = (dynamic_cast<TDirectory*>(obj))->GetList();
TIter it(listSub);
while (TObject* obj1 = it())
writedir2current(obj1);
cur->cd();
}
}
// ---------------------------------------------------------------------------------------------------------------------
//
void CbmL1::GenerateMaterialMaps()
{
LOG(info) << "Generating material maps...";
auto timerStart = std::chrono::high_resolution_clock::now();
cbm::ca::tools::MaterialHelper matHelper;
matHelper.SetSafeMaterialInitialization(fDoSafeMaterialInitialization);
if (!fMatBudgetParallelProjection) {
matHelper.SetDoRadialProjection(fTargetZ);
}
matHelper.SetNraysPerDim(fMatBudgetNrays);
double zLast = fTargetZ + 1.; // some gap (+1cm) to skip the target material
std::vector<ca::StationInitializer*> vpActiveStations;
vpActiveStations.reserve(fInitManager.GetNstationsActive());
for (auto& station : fInitManager.GetStationInfo()) { // loop over active + inactive station
if (station.GetTrackingStatus()) {
vpActiveStations.push_back(&station); }
}
LOG(info) << "Generating material maps for stations: ";
for (const auto* pSta : vpActiveStations) {
LOG(info) << "\t- z = " << pSta->GetZref() << " cm";
}
for (unsigned int ist = 0; ist < vpActiveStations.size(); ++ist) {
auto* pStation = vpActiveStations[ist];
double z1 = pStation->GetZmax();
double z2 = z1;
if (ist < vpActiveStations.size() - 1) {
// split materials between the stations at the middle
auto* pStationNext = vpActiveStations[ist + 1];
z2 = pStationNext->GetZmin();
}
double zNew = 0.5 * (z1 + z2);
double maxXY = 1.3 * std::max(std::fabs(pStation->GetXmax()), std::fabs(pStation->GetYmax()));
//double maxXY = 80;
// calculate n bins from the minimal pitch
int nBins = static_cast<int>(std::ceil(2. * maxXY / fMatBudgetPitch));
if (nBins < 1) {
LOG(fatal) << " material nBins " << nBins << " is not positive, something is wrong";
}
if (nBins > fMatBudgetNbins) {
nBins = fMatBudgetNbins;
}
ca::MaterialMap matBudget = matHelper.GenerateMaterialMap(pStation->GetZref(), zLast, zNew, maxXY, nBins);
pStation->SetMaterialMap(matBudget);
LOG(info) << "Generated material map for tracking station " << ist << " at z = " << pStation->GetZref() << " cm."
<< " Material is collected between z = " << zLast << " and z = " << zNew;
zLast = zNew;
}
auto timerEnd = std::chrono::high_resolution_clock::now();
double materialGenerationTime = (double) (std::chrono::duration<double>(timerEnd - timerStart).count());
LOG(info) << "Generating material maps... done! (it took " << materialGenerationTime << " s)";
}
// ---------------------------------------------------------------------------------------------------------------------
//
void CbmL1::IdealTrackFinder()
{
fpAlgo->fRecoTracks.clear();
fpAlgo->fRecoHits.clear();
for (auto& MC : fMCData.GetTrackContainer()) {
if (!MC.IsReconstructable()) continue;
if (!(MC.GetId() >= 0)) continue;
if (MC.GetNofHits() < 4) continue;
CaTrack algoTr;
algoTr.fNofHits = 0;
ca::Vector<int> hitIndices("CbmL1::hitIndices", fpAlgo->GetParameters().GetNstationsActive(), -1);
for (unsigned int iH : MC.GetHitIndexes()) {
const CbmL1HitDebugInfo& hit = fvHitDebugInfo[iH];
const int iStation = hit.GetStationId();
if (iStation >= 0) hitIndices[iStation] = iH;
}
for (int iH = 0; iH < fpAlgo->GetParameters().GetNstationsActive(); iH++) {
const int hitI = hitIndices[iH];
if (hitI < 0) continue;
// fpAlgo->fRecoHits.push_back(hitI);
algoTr.fNofHits++;
}
if (algoTr.fNofHits < 3) continue;
for (int iH = 0; iH < fpAlgo->GetParameters().GetNstationsActive(); iH++) {
const int hitI = hitIndices[iH];
if (hitI < 0) continue;
fpAlgo->fRecoHits.push_back(hitI);
}
algoTr.fParFirst.X() = MC.GetStartX();
algoTr.fParFirst.Y() = MC.GetStartY();
algoTr.fParFirst.Z() = MC.GetStartZ();
algoTr.fParFirst.Tx() = MC.GetTx();
algoTr.fParFirst.Ty() = MC.GetTy();
algoTr.fParFirst.Qp() = MC.GetCharge() / MC.GetP();
algoTr.fParFirst.Time() = MC.GetStartT();
fpAlgo->fRecoTracks.push_back(algoTr);
}
} // void CbmL1::IdealTrackFinder()
/// ----- STandAlone Package service-functions -----------------------------
// ---------------------------------------------------------------------------------------------------------------------
//
void CbmL1::DefineSTAPNames(const char* dirName)
{
// FIXME: SZh 01.03.2023: Clean STAP names
namespace bfs = boost::filesystem;
if (fSTAPDataMode == 0) {
return;
}
// Define file prefix (/path/to/data/setup.reco.root -> "setup.reco")
bfs::path pathToRecoOutput = FairRunAna::Instance()->GetUserOutputFileName().Data();
fSTAPDataPrefix = pathToRecoOutput.filename().string();
fSTAPDataPrefix.ReplaceAll(".root", "");
fSTAPDataPrefix = fSTAPDataPrefix.Strip(TString::EStripType::kBoth, '.');
TString sDirName = TString(dirName);
if (sDirName.Length() == 0) {
fSTAPDataDir = pathToRecoOutput.parent_path().string();
}
else if (bfs::exists(sDirName.Data()) && bfs::is_directory(sDirName.Data())) {
fSTAPDataDir = sDirName;
}
else {
fSTAPDataDir = ".";
}
// TODO: Clean-up the names and pathes
if (fSTAPDataDir.Length() == 0) {
fSTAPDataDir = ".";
}
LOG(info) << "CbmL1: STAP data root directory is \033[1;32m" << bfs::system_complete(fSTAPDataDir.Data())
<< "\033[0m";
if (fSTAPDataMode == 4) {
return;
}
// Directory for handling L1InputData objects
TString sInputDataDir = fSTAPDataDir + "/" + fSTAPDataPrefix + "_cainputdata";
if (!bfs::exists(sInputDataDir.Data())) {
LOG(warn) << "CbmL1: directory for tracking input data does not exist. It will be created";
bfs::create_directories(sInputDataDir.Data());
}
LOG(info) << "CbmL1: STAP tracking input jobs directory is \033[1;32m" << bfs::system_complete(sInputDataDir.Data())
<< "\033[0m";
}
// ---------------------------------------------------------------------------------------------------------------------
//
void CbmL1::WriteSTAPParamObject()
{
TString filename = fSTAPParamFile.IsNull()
? fSTAPDataDir + "/" + fSTAPDataPrefix + "." + TString(kSTAPParamSuffix.data())
: fSTAPParamFile;
fInitManager.WriteParametersObject(filename.Data());
}
// ---------------------------------------------------------------------------------------------------------------------
//
void CbmL1::WriteSTAPAlgoInputData(int iJob) // must be called after ReadEvent
{
TString filename =
fSTAPDataDir + "/" + fSTAPDataPrefix + "_cainputdata/" + +TString::Format(kSTAPAlgoIDataSuffix.data(), iJob);
// Write file
fpIODataManager->WriteInputData(filename.Data());
}
// ---------------------------------------------------------------------------------------------------------------------
//
void CbmL1::WriteSTAPPerfInputData() // must be called after ReadEvent
{
LOG(fatal) << "CbmL1: Running in standalone mode is not available at the moment. It will be updated soon...";
}
// ---------------------------------------------------------------------------------------------------------------------
//
void CbmL1::ReadSTAPParamObject()
{
TString filename = fSTAPDataDir + "/" + fSTAPDataPrefix + "." + TString(kSTAPParamSuffix.data());
fInitManager.ReadParametersObject(filename.Data());
}
// ---------------------------------------------------------------------------------------------------------------------
//
void CbmL1::ReadSTAPAlgoInputData(int iJob)
{
TString filename = fSTAPDataDir + "/" + TString(kSTAPAlgoIDataDir) + "/" + fSTAPDataPrefix + "."
+ TString::Format(kSTAPAlgoIDataSuffix.data(), iJob);
// Read file
fpIODataManager->ReadInputData(filename.Data());
}
// ---------------------------------------------------------------------------------------------------------------------
//
void CbmL1::ReadSTAPPerfInputData()
{
LOG(fatal) << "CbmL1: Running in standalone mode is not available at the moment. It will be updated soon...";
}
// ---------------------------------------------------------------------------------------------------------------------
//
void CbmL1::DumpMaterialToFile(TString fileName)
{
TFile* f = new TFile(fileName, "RECREATE"); f->cd();
const auto& param = fpAlgo->GetParameters();
const auto& map = param.GetThicknessMaps();
for (int ist = 0; ist < param.GetNstationsActive(); ist++) {
//const ca::Station& st = param.GetStation(ist);
const auto& mat = map[ist];
TString name = Form("tracking_station%d", ist);
TString title = Form("Tracking station %d: Rad. thickness in [%s]. Z region [%.2f,%.2f] cm.", ist, "%",
mat.GetZmin(), mat.GetZmax());
if (fMatBudgetParallelProjection) {
title += " Horisontal projection.";
}
else {
title += " Radial projection.";
}
TH2F* h = new TH2F(name.Data(), title.Data(), mat.GetNbins(), -mat.GetXYmax(), mat.GetXYmax(), mat.GetNbins(),
-mat.GetXYmax(), mat.GetXYmax());
h->GetXaxis()->SetTitle("X [cm]");
h->GetYaxis()->SetTitle("Y [cm]");
for (int ix = 0; ix < mat.GetNbins(); ix++) {
for (int iy = 0; iy < mat.GetNbins(); iy++) {
h->SetBinContent(ix + 1, iy + 1, 100. * mat.GetRadThickBin(ix, iy));
}
}
}
f->Write();
}
// ---------------------------------------------------------------------------------------------------------------------
//
void CbmL1::GetTargetInfo()
{
// Loop over all nodes till a node with name "target" is found
// Extract the required infrmation from the node and store it in the
// proper structure
// The complete logic depends on the naming convention of the target.
// If the node doesn't contain the string target the procedure will fail
fTargetX = 1.e10;
fTargetY = 1.e10;
fTargetZ = 1.e10;
TString targetPath;
TGeoNode* targetNode{nullptr};
FindTargetNode(targetPath, targetNode);
if (!targetNode) {
LOG(fatal) << "L1: can not find the target!";
}
Double_t local[3] = {0., 0., 0.}; // target centre, local c.s.
Double_t global[3]; // target centre, global c.s.
gGeoManager->cd(targetPath);
gGeoManager->GetCurrentMatrix()->LocalToMaster(local, global);
fTargetX = global[0];
fTargetY = global[1];
fTargetZ = global[2];
LOG(info) << " Target found: \"" << targetPath << "\" at ( " << fTargetX << " " << fTargetY << " " << fTargetZ
<< " ) ";
}
// ---------------------------------------------------------------------------------------------------------------------
// Target finder - recursive routine
//
void CbmL1::FindTargetNode(TString& targetPath, TGeoNode*& targetNode)
{
if (!targetNode) { // init at the top of the tree
targetNode = gGeoManager->GetTopNode();
targetPath = "/" + TString(targetNode->GetName());
}
if (TString(targetNode->GetName()).Contains("target")) {
return;
}
for (Int_t iNode = 0; iNode < targetNode->GetNdaughters(); iNode++) {
TGeoNode* newNode = targetNode->GetDaughter(iNode);
TString newPath = targetPath + "/" + newNode->GetName();
FindTargetNode(newPath, newNode);
if (newNode) {
targetPath = newPath;
targetNode = newNode;
return;
}
}
targetPath = "";
targetNode = nullptr;
}