From 2e7740669cbaed9fa4bb88af184b6bd6c6144870 Mon Sep 17 00:00:00 2001
From: "s.zharko@gsi.de" <s.zharko@gsi.de>
Date: Tue, 5 Dec 2023 20:51:49 +0100
Subject: [PATCH] CA: CbmL1 clean-up
---
algo/ca/core/pars/CaMaterialMap.cxx | 11 +
algo/ca/core/pars/CaMaterialMap.h | 3 +
algo/ca/core/pars/CaParameters.cxx | 28 +-
algo/ca/core/pars/CaStationInitializer.cxx | 5 +
reco/L1/CbmL1.cxx | 427 ++++--------------
reco/L1/CbmL1.h | 103 +----
reco/L1/CbmL1MCTrack.cxx | 21 +-
reco/L1/CbmL1Performance.cxx | 485 ---------------------
reco/L1/L1Algo/utils/L1AlgoPulls.cxx | 2 +-
reco/L1/L1Algo/utils/L1AlgoPulls.h | 7 +
reco/L1/qa/CbmCaOutputQa.cxx | 54 ---
11 files changed, 137 insertions(+), 1009 deletions(-)
diff --git a/algo/ca/core/pars/CaMaterialMap.cxx b/algo/ca/core/pars/CaMaterialMap.cxx
index 8b02c09dd4..75fc25c572 100644
--- a/algo/ca/core/pars/CaMaterialMap.cxx
+++ b/algo/ca/core/pars/CaMaterialMap.cxx
@@ -133,3 +133,14 @@ void MaterialMap::Swap(MaterialMap& other) noexcept
std::swap(fZmax, other.fZmax);
std::swap(fTable, other.fTable); // Probably can cause segmentation violation (did not understand)
}
+
+//------------------------------------------------------------------------------------------------------------------------------------
+//
+std::string MaterialMap::ToString() const
+{
+ std::stringstream msg;
+ using std::setw;
+ msg << "[mat.map] nBins " << setw(6) << fNbins << ", XYmax " << setw(12) << fXYmax << " , z (ref, min, max) "
+ << setw(12) << fZref << ' ' << setw(12) << fZmin << ' ' << setw(12) << fZmax;
+ return msg.str();
+}
diff --git a/algo/ca/core/pars/CaMaterialMap.h b/algo/ca/core/pars/CaMaterialMap.h
index 9d80f2d6c3..9b95e7b319 100644
--- a/algo/ca/core/pars/CaMaterialMap.h
+++ b/algo/ca/core/pars/CaMaterialMap.h
@@ -103,6 +103,9 @@ namespace cbm::algo::ca
/// \brief Get bin index for (x,y). Returns -1 when outside of the map
int GetBin(float x, float y) const;
+ /// \brief String representation of the object
+ std::string ToString() const;
+
private:
int fNbins = constants::Undef<int>; ///< Number of rows (== N columns) in the material budget table
float fXYmax = constants::Undef<float>; ///< Size of the station in x and y dimensions [cm]
diff --git a/algo/ca/core/pars/CaParameters.cxx b/algo/ca/core/pars/CaParameters.cxx
index 516e86f3d7..a6658f206d 100644
--- a/algo/ca/core/pars/CaParameters.cxx
+++ b/algo/ca/core/pars/CaParameters.cxx
@@ -255,15 +255,15 @@ std::string Parameters::ToString(int verbosity, int indentLevel) const
msg << indent << indentCh << clrs::CLb << "NUMBER OF STATIONS:\n" << clrs::CL;
msg << indent << indentCh << indentCh << "Number of stations (Geometry): ";
for (int iDet = 0; iDet < constants::size::MaxNdetectors; ++iDet) {
- msg << setw(2) << setfill(' ') << this->GetNstationsGeometry(static_cast<EDetectorID>(iDet)) << ' ';
+ msg << setw(2) << this->GetNstationsGeometry(static_cast<EDetectorID>(iDet)) << ' ';
}
- msg << " | total = " << setw(2) << setfill(' ') << this->GetNstationsGeometry();
+ msg << " | total = " << setw(2) << this->GetNstationsGeometry();
msg << '\n';
msg << indent << indentCh << indentCh << "Number of stations (Active): ";
for (int iDet = 0; iDet < constants::size::MaxNdetectors; ++iDet) {
msg << setw(2) << setfill(' ') << this->GetNstationsActive(static_cast<EDetectorID>(iDet)) << ' ';
}
- msg << " | total = " << setw(2) << setfill(' ') << this->GetNstationsActive();
+ msg << " | total = " << setw(2) << this->GetNstationsActive();
msg << '\n' << indent << indentCh << indentCh << clrs::CL << "Geometry station indices: ";
for (int iStGeo = 0; iStGeo < this->GetNstationsGeometry(); ++iStGeo) {
bool isActive = fvGeoToActiveMap[iStGeo] != -1;
@@ -287,13 +287,31 @@ std::string Parameters::ToString(int verbosity, int indentLevel) const
msg << clrs::CL << '\n';
msg << indent << indentCh << clrs::CLb << "STATIONS:\n" << clrs::CL;
- msg << indent << indentCh << setw(9) << setfill(' ') << "Active ID" << ' ';
+ msg << indent << indentCh << setw(9) << "Active ID" << ' ';
msg << fStations[0].ToString(1, 0, true) << '\n';
for (int iStAct = 0; iStAct < this->GetNstationsActive(); ++iStAct) {
- msg << indent << indentCh << setw(9) << setfill(' ') << iStAct << ' ';
+ msg << indent << indentCh << setw(9) << iStAct << ' ';
msg << fStations[iStAct].ToString(verbosity, 0) << '\n';
}
+ msg << indent << indentCh << clrs::CLb << "STATION MATERIAL MAPS:\n" << clrs::CL;
+ msg << indent << indentCh;
+ msg << setw(9) << "Active ID" << ' ';
+ msg << setw(9) << "N bins" << ' ';
+ msg << setw(11) << "xy max.[cm]" << ' ';
+ msg << setw(10) << "z min.[cm]" << ' ';
+ msg << setw(10) << "z max.[cm]" << ' ';
+ msg << setw(10) << "z ref.[cm]" << '\n';
+ for (int iStAct = 0; iStAct < this->GetNstationsActive(); ++iStAct) {
+ const auto& matMap = fThickMap[iStAct];
+ msg << indent << indentCh << setw(9) << iStAct << ' ';
+ msg << setw(9) << matMap.GetNbins() << ' ';
+ msg << setw(11) << matMap.GetXYmax() << ' ';
+ msg << setw(10) << matMap.GetZmin() << ' ';
+ msg << setw(10) << matMap.GetZmax() << ' ';
+ msg << setw(10) << matMap.GetZref() << '\n';
+ }
+
msg << indent << clrs::CLb << "DEV FLAGS:" << clrs::CL << " (for debug only)\n";
msg << indent << indentCh << "Hits search area is ignored: " << fDevIsIgnoreHitSearchAreas << '\n';
msg << indent << indentCh << "Non-approx. field is used: " << fDevIsMatchDoubletsViaMc << '\n';
diff --git a/algo/ca/core/pars/CaStationInitializer.cxx b/algo/ca/core/pars/CaStationInitializer.cxx
index b6679f9665..691fa04968 100644
--- a/algo/ca/core/pars/CaStationInitializer.cxx
+++ b/algo/ca/core/pars/CaStationInitializer.cxx
@@ -273,6 +273,11 @@ void StationInitializer::SetTrackingStatus(bool flag)
{
fTrackingStatus = flag;
fInitController.SetFlag(EInitKey::kTrackingStatus);
+ if (!fTrackingStatus && !fInitController.GetFlag(EInitKey::kThicknessMap)) {
+ // Set initialized status for the material map, if the station is not active, since its material is accounted in the
+ // previous tracking station.
+ fInitController.SetFlag(EInitKey::kThicknessMap);
+ }
}
// ---------------------------------------------------------------------------------------------------------------------
diff --git a/reco/L1/CbmL1.cxx b/reco/L1/CbmL1.cxx
index b1b8ab0048..8b1e545a7e 100644
--- a/reco/L1/CbmL1.cxx
+++ b/reco/L1/CbmL1.cxx
@@ -41,7 +41,6 @@
#include "CbmStsFindTracks.h"
#include "CbmStsHit.h"
#include "CbmTrackingDetectorInterfaceInit.h"
-#include "FairEventHeader.h"
#include "FairField.h"
#include "FairRunAna.h"
#include "TGeoArb8.h"
@@ -198,25 +197,23 @@ InitStatus CbmL1::Init()
fUseTRD = false;
fUseTOF = false;
- FairRootManager* fairManager = FairRootManager::Instance();
-
if (ca::Framework::TrackingMode::kSts == fTrackingMode) {
- fUseMVD = 1;
- fUseSTS = 1;
- fUseMUCH = 0;
- fUseTRD = 0;
- fUseTOF = 0;
+ 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::Framework::TrackingMode::kMcbm == fTrackingMode) {
- fUseMVD = 1;
- fUseSTS = 1;
- fUseMUCH = 1;
- fUseTRD = 1;
- fUseTOF = 1;
+ fUseMVD = true;
+ fUseSTS = true;
+ fUseMUCH = true;
+ fUseTRD = true;
+ fUseTOF = true;
TString tag;
CbmSetup::Instance()->GetGeoTag(ECbmModuleId::kMuch, tag);
if (tag.Contains("mcbm")) { // currently disable tracking in much for all mcbm setups
@@ -227,11 +224,11 @@ InitStatus CbmL1::Init()
if (ca::Framework::TrackingMode::kGlobal == fTrackingMode) {
//at the moment trd2d tracking only
- fUseMVD = 1;
- fUseSTS = 1;
- fUseMUCH = 1;
- fUseTRD = 1;
- fUseTOF = 0;
+ fUseMVD = true;
+ fUseSTS = true;
+ fUseMUCH = true;
+ fUseTRD = true;
+ fUseTOF = false;
fInitManager.DevSetUseOfOriginalField();
fInitManager.DevSetIgnoreHitSearchAreas(true);
@@ -247,167 +244,6 @@ InitStatus CbmL1::Init()
CheckDetectorPresence();
- fpStsPoints = nullptr;
- fpMvdPoints = nullptr;
- fpMuchPoints = nullptr;
- fpTrdPoints = nullptr;
- fpTofPoints = nullptr;
- fpMcTracks = nullptr;
-
- fpMvdHitMatches = nullptr;
- fpTrdHitMatches = nullptr;
- fpMuchHitMatches = nullptr;
- fpTofHitMatches = nullptr;
-
- fpStsClusters = nullptr;
-
- fTimeSlice = (CbmTimeSlice*) fairManager->GetObject("TimeSlice.");
- if (!fTimeSlice) {
- LOG(fatal) << GetName() << ": No time slice branch in the tree!";
- }
-
- fpStsClusters = dynamic_cast<TClonesArray*>(fairManager->GetObject("StsCluster"));
-
- fpStsHits = dynamic_cast<TClonesArray*>(fairManager->GetObject("StsHit"));
-
- if (!fUseMUCH) {
- fpMuchPixelHits = nullptr;
- fpMuchDigis = nullptr; // NOTE: Used only with fPerformance = true
- fpMuchDigiMatches = nullptr; // NOTE: Used only with fPerformance = true
- fpMuchClusters = nullptr; // NOTE: Used only with fPerformance = true
- fpMuchPoints = nullptr; // NOTE: Used only with fPerformance = true
- fpMuchHitMatches = nullptr; // NOTE: Used only with fPerformance = true
- }
- else {
- fpMuchPixelHits = (TClonesArray*) fairManager->GetObject("MuchPixelHit");
- }
-
- if (!fUseTRD) {
- fpTrdPoints = nullptr;
- fpTrdHitMatches = nullptr;
- fpTrdHits = nullptr;
- }
- else {
- fpTrdHits = dynamic_cast<TClonesArray*>(fairManager->GetObject("TrdHit"));
- }
-
- if (!fUseTOF) {
- fpTofPoints = nullptr;
- fpTofHitMatches = nullptr;
- fpTofHits = nullptr;
- }
- else {
- fpTofHits = (TClonesArray*) fairManager->GetObject("TofHit");
- }
-
- // misalignment
- for (int i = 0; i != static_cast<int>(cbm::algo::ca::EDetectorID::kEND); ++i) {
- auto detId = cbm::algo::ca::EDetectorID(i);
- fInitManager.SetMisalignmentX(detId, fvMisalignment[detId][0]);
- fInitManager.SetMisalignmentY(detId, fvMisalignment[detId][1]);
- fInitManager.SetMisalignmentT(detId, fvMisalignment[detId][2]);
- LOG(info) << "CbmL1: misalignment for " << cbm::ca::kDetName[detId] << " is set to " << fvMisalignment[detId][0]
- << " " << fvMisalignment[detId][1] << " " << fvMisalignment[detId][2];
- }
-
- if (fPerformance) {
-
- CbmMCDataManager* mcManager = (CbmMCDataManager*) fairManager->GetObject("MCDataManager");
- if (nullptr == mcManager) LOG(fatal) << GetName() << ": No CbmMCDataManager!";
-
- fpMcEventHeader = mcManager->GetObject("MCEventHeader.");
-
- fpMcTracks = mcManager->InitBranch("MCTrack");
-
- if (nullptr == fpMcTracks) LOG(fatal) << GetName() << ": No MCTrack data!";
- if (nullptr == fpMcEventHeader) LOG(fatal) << GetName() << ": No MC event header data!";
-
- fpMcEventList = (CbmMCEventList*) fairManager->GetObject("MCEventList.");
- if (nullptr == fpMcEventList) LOG(fatal) << GetName() << ": No MCEventList data!";
-
- if (fUseMVD) {
- fpMvdPoints = mcManager->InitBranch("MvdPoint");
- fpMvdDigiMatches = dynamic_cast<TClonesArray*>(fairManager->GetObject("MvdDigiMatch"));
- fpMvdHitMatches = dynamic_cast<TClonesArray*>(fairManager->GetObject("MvdHitMatch"));
- if (!fpMvdHitMatches) {
- LOG(error) << "No fpMvdHitMatches provided, performance is not done correctly";
- }
- }
-
- if (fUseSTS) {
- fpStsPoints = mcManager->InitBranch("StsPoint");
- fpStsHitMatches = dynamic_cast<TClonesArray*>(fairManager->GetObject("StsHitMatch"));
- fpStsClusterMatches = dynamic_cast<TClonesArray*>(fairManager->GetObject("StsClusterMatch"));
- if (nullptr == fpStsPoints) LOG(fatal) << GetName() << ": No StsPoint data!";
- }
-
- fpTrdPoints = mcManager->InitBranch("TrdPoint");
-
- if (!fUseTRD) {
- fpTrdHitMatches = nullptr;
- }
- else {
- fpTrdHitMatches = (TClonesArray*) fairManager->GetObject("TrdHitMatch");
- }
-
- if (!fUseMUCH) {
- fpMuchPoints = nullptr;
- fpMuchHitMatches = nullptr;
- }
- else {
- fpMuchDigis = (TClonesArray*) fairManager->GetObject("MuchDigi");
- fpMuchDigiMatches = (TClonesArray*) fairManager->GetObject("MuchDigiMatch");
- fpMuchClusters = (TClonesArray*) fairManager->GetObject("MuchCluster");
- fpMuchPoints = mcManager->InitBranch("MuchPoint");
- fpMuchHitMatches = dynamic_cast<TClonesArray*>(fairManager->GetObject("MuchPixelHitMatch"));
- }
-
- fpTofPoints = mcManager->InitBranch("TofPoint");
-
- if (!fUseTOF) {
- fpTofHitMatches = nullptr;
- }
- else {
- fpTofHitMatches = static_cast<TClonesArray*>(fairManager->GetObject("TofHitMatch"));
- }
-
- // Check, if performance is doable
- auto CheckBranch = [&](const std::string& name, const auto* ptr) {
- if (!ptr) {
- LOG(error) << GetName() << "Branch " << name << " not found, performance is set 0";
- fPerformance = 0;
- }
- };
-
- if (fUseMVD) {
- CheckBranch("MvdPoint", fpMvdPoints);
- CheckBranch("MvdHitMatches", fpMvdHitMatches);
- }
- if (fUseSTS) {
- CheckBranch("StsPoint", fpStsPoints);
- CheckBranch("StsHitMatches", fpStsHitMatches);
- }
- if (fUseMUCH) {
- CheckBranch("MuchPoint", fpMuchPoints);
- CheckBranch("MuchPixelHitMatch", fpMuchHitMatches);
- }
- if (fUseTRD) {
- CheckBranch("TrdPoint", fpTrdPoints);
- CheckBranch("TrdHitMatch", fpTrdHitMatches);
- }
- if (fUseTOF) {
- CheckBranch("TofPoint", fpTofPoints);
- CheckBranch("TofHitMatch", fpTofHitMatches);
- }
- }
-
- if (!fUseMVD) {
- fpMvdHits = nullptr;
- }
- else {
- fpMvdHits = dynamic_cast<TClonesArray*>(fairManager->GetObject("MvdHit"));
- }
-
// *****************************
// ** **
// ** GEOMETRY INITIALIZATION **
@@ -417,13 +253,6 @@ InitStatus CbmL1::Init()
// Read parameters object from a binary file
if (2 == fSTAPDataMode) {
this->ReadSTAPParamObject();
-
- fNMvdStationsGeom = fInitManager.GetNstationsGeometry(ca::EDetectorID::kMvd);
- fNStsStationsGeom = fInitManager.GetNstationsGeometry(ca::EDetectorID::kSts);
- fNTrdStationsGeom = fInitManager.GetNstationsGeometry(ca::EDetectorID::kTrd);
- fNMuchStationsGeom = fInitManager.GetNstationsGeometry(ca::EDetectorID::kMuch);
- fNTofStationsGeom = fInitManager.GetNstationsGeometry(ca::EDetectorID::kTof);
- fNStationsGeom = fInitManager.GetNstationsGeometry();
}
// Parameters initialization, based on the FairRunAna chain
else {
@@ -462,12 +291,11 @@ InitStatus CbmL1::Init()
}
}
- fNMvdStationsGeom = (fUseMVD) ? mvdInterface->GetNtrackingStations() : 0;
- fNStsStationsGeom = (fUseSTS) ? stsInterface->GetNtrackingStations() : 0;
- fNMuchStationsGeom = (fUseMUCH) ? muchInterface->GetNtrackingStations() : 0;
- fNTrdStationsGeom = (fUseTRD) ? trdInterface->GetNtrackingStations() : 0;
- fNTofStationsGeom = (fUseTOF) ? tofInterface->GetNtrackingStations() : 0;
- fNStationsGeom = fNMvdStationsGeom + fNStsStationsGeom + fNMuchStationsGeom + fNTrdStationsGeom + fNTofStationsGeom;
+ 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 **
@@ -502,6 +330,7 @@ InitStatus CbmL1::Init()
// 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 **
@@ -532,6 +361,19 @@ InitStatus CbmL1::Init()
}
fInitManager.SetActiveDetectorIDs(vActiveTrackingDetectorIDs);
+ // *********************************
+ // ** Misalignment initialization **
+ // *********************************
+
+ for (int iDet = 0; iDet != static_cast<int>(cbm::algo::ca::EDetectorID::kEND); ++iDet) {
+ auto detId = cbm::algo::ca::EDetectorID(iDet);
+ fInitManager.SetMisalignmentX(detId, fvMisalignment[detId][0]);
+ fInitManager.SetMisalignmentY(detId, fvMisalignment[detId][1]);
+ fInitManager.SetMisalignmentT(detId, fvMisalignment[detId][2]);
+ LOG(info) << "CbmL1: misalignment for " << cbm::ca::kDetName[detId] << " is set to " << fvMisalignment[detId][0]
+ << " " << fvMisalignment[detId][1] << " " << fvMisalignment[detId][2];
+ }
+
// *************************************
// ** Stations layout initialization **
// *************************************
@@ -542,7 +384,7 @@ InitStatus CbmL1::Init()
// *** MVD stations info ***
if (fUseMVD) {
- for (int iSt = 0; iSt < fNMvdStationsGeom; ++iSt) {
+ 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
@@ -565,7 +407,7 @@ InitStatus CbmL1::Init()
// *** STS stations info ***
if (fUseSTS) {
- for (int iSt = 0; iSt < fNStsStationsGeom; ++iSt) {
+ 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
@@ -589,7 +431,7 @@ InitStatus CbmL1::Init()
// *** MuCh stations info ***
if (fUseMUCH) {
- for (int iSt = 0; iSt < fNMuchStationsGeom; ++iSt) {
+ 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
@@ -613,7 +455,7 @@ InitStatus CbmL1::Init()
// *** TRD stations info ***
if (fUseTRD) {
- for (int iSt = 0; iSt < fNTrdStationsGeom; ++iSt) {
+ 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((iSt == 1 || iSt == 3) ? 6 : 3); // MuCh
@@ -640,7 +482,7 @@ InitStatus CbmL1::Init()
// *** TOF stations info ***
if (fUseTOF) {
- for (int iSt = 0; iSt < fNTofStationsGeom; ++iSt) {
+ 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);
@@ -749,20 +591,14 @@ InitStatus CbmL1::Init()
LOG(info) << fpParameters->ToString(1);
LOG(info) << "----- Numbers of stations active in tracking -----";
- fNMvdStations = fpParameters->GetNstationsActive(ca::EDetectorID::kMvd);
- fNStsStations = fpParameters->GetNstationsActive(ca::EDetectorID::kSts);
- fNTrdStations = fpParameters->GetNstationsActive(ca::EDetectorID::kTrd);
- fNMuchStations = fpParameters->GetNstationsActive(ca::EDetectorID::kMuch);
- fNTofStations = fpParameters->GetNstationsActive(ca::EDetectorID::kTof);
- fNStations = fpParameters->GetNstationsActive();
-
+ LOG(info) << " MVD: " << fpParameters->GetNstationsActive(ca::EDetectorID::kMvd);
+ LOG(info) << " STS: " << fpParameters->GetNstationsActive(ca::EDetectorID::kSts);
+ LOG(info) << " MuCh: " << fpParameters->GetNstationsActive(ca::EDetectorID::kTrd);
+ LOG(info) << " TRD: " << fpParameters->GetNstationsActive(ca::EDetectorID::kMuch);
+ LOG(info) << " TOF: " << fpParameters->GetNstationsActive(ca::EDetectorID::kTof);
+ LOG(info) << " Total: " << fpParameters->GetNstationsActive();
- LOG(info) << " MVD: " << fNMvdStations;
- LOG(info) << " STS: " << fNStsStations;
- LOG(info) << " MuCh: " << fNMuchStations;
- LOG(info) << " TRD: " << fNTrdStations;
- LOG(info) << " ToF: " << fNTofStations;
- LOG(info) << " Total: " << fNStations;
+ fNStations = fpParameters->GetNstationsActive();
// monitor the material
{
@@ -787,28 +623,19 @@ InitStatus CbmL1::Init()
void CbmL1::Reconstruct(CbmEvent* event)
{
- fpTSReader->ReadEvent(event);
- LOG(info) << "READ: hit ids = " << fvExternalHits.size() << ", hits = " << fpIODataManager->GetNofHits()
- << ", dbg hits = " << fvHitDebugInfo.size();
-
-
ca::TrackingMonitorData monitorData;
- LOG_IF(info, fVerbose > 1) << "\nCbmL1::Exec " << (event ? "event " : "time slice ") << fEventNo << " ...\n";
-
- // ----- Read data from branches and send data from IODataManager to ca::Framework ----------------------------------------
-
-
+ 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 (ca::HitIndex_t i = 0; i < fpAlgo->GetInputData().GetNhits(); i++) {
- const ca::Hit& h = fpAlgo->GetInputData().GetHit(i);
- fMaterialMonitor[h.Station()].MarkActiveBin(h.X(), h.Y());
+ 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()
@@ -816,7 +643,6 @@ void CbmL1::Reconstruct(CbmEvent* event)
// FieldApproxCheck();
// FieldIntegralCheck();
- fpAlgo->ReceiveInputData(fpIODataManager->TakeInputData());
fpAlgo->SetMonitorData(monitorData);
LOG_IF(info, fVerbose > 1) << "L1 Track finder...";
@@ -836,19 +662,20 @@ void CbmL1::Reconstruct(CbmEvent* event)
int trackFirstHit = 0;
// FIXME: Rewrite
- for (auto it = fpAlgo->fRecoTracks.begin(); it != fpAlgo->fRecoTracks.end(); trackFirstHit += it->fNofHits, it++) {
+ for (const auto& caTrk : fpAlgo->fRecoTracks) {
CbmL1Track t;
- t.Set(it->fParFirst);
- t.TLast.Set(it->fParLast);
- t.Tpv.Set(it->fParPV);
+ t.Set(caTrk.fParFirst);
+ t.TLast.Set(caTrk.fParLast);
+ t.Tpv.Set(caTrk.fParPV);
t.Hits.clear();
- for (int i = 0; i < it->fNofHits; i++) {
+ 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);
}
@@ -857,9 +684,7 @@ void CbmL1::Reconstruct(CbmEvent* event)
// output performance
if (fPerformance) {
- if (fVerbose > 1) {
- LOG(info) << "Performance...";
- }
+ LOG_IF(info, fVerbose) << "Performance...";
fpMCModule->MatchTracks(); // matches reco and MC tracks, fills the MC-truth fields of reco tracks
@@ -982,19 +807,26 @@ void CbmL1::GenerateMaterialMaps()
double zLast = fTargetZ + 1.; // some gap (+1cm) to skip the target material
- auto& vStations = fInitManager.GetStationInfo();
- for (unsigned int ist = 0; ist < vStations.size(); ist++) {
- auto& station = vStations[ist];
- double z1 = station.GetZmax();
+ std::vector<ca::StationInitializer*> vpActiveStations;
+ vpActiveStations.reserve(fInitManager.GetNstationsActive());
+ for (auto& station : fInitManager.GetStationInfo()) { // loop over active + inactive stations
+ if (station.GetTrackingStatus()) {
+ vpActiveStations.push_back(&station);
+ }
+ }
+
+ for (unsigned int ist = 0; ist < vpActiveStations.size(); ++ist) {
+ auto* pStation = vpActiveStations[ist];
+ double z1 = pStation->GetZmax();
double z2 = z1;
- if (ist < vStations.size() - 1) {
+ if (ist < vpActiveStations.size() - 1) {
// split materials between the stations at the middle
- auto& stationNext = vStations[ist + 1];
- z2 = stationNext.GetZmin();
+ auto* pStationNext = vpActiveStations[ist + 1];
+ z2 = pStationNext->GetZmin();
}
double zNew = 0.5 * (z1 + z2);
- double maxXY = 1.3 * std::max(station.GetXmax(), station.GetYmax());
+ double maxXY = 1.3 * std::max(pStation->GetXmax(), pStation->GetYmax());
//double maxXY = 80;
// calculate n bins from the minimal pitch
int nBins = static_cast<int>(std::ceil(2. * maxXY / fMatBudgetPitch));
@@ -1005,11 +837,11 @@ void CbmL1::GenerateMaterialMaps()
nBins = fMatBudgetNbins;
}
- ca::MaterialMap matBudget = matHelper.GenerateMaterialMap(station.GetZref(), zLast, zNew, maxXY, nBins);
+ ca::MaterialMap matBudget = matHelper.GenerateMaterialMap(pStation->GetZref(), zLast, zNew, maxXY, nBins);
- station.SetMaterialMap(matBudget);
+ pStation->SetMaterialMap(matBudget);
- LOG(info) << "Generated material map for tracking station " << ist << " at z = " << station.GetZref() << " cm."
+ 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;
@@ -1175,113 +1007,6 @@ void CbmL1::ReadSTAPPerfInputData()
// ---------------------------------------------------------------------------------------------------------------------
//
-void CbmL1::WriteSIMDKFData()
-{
- // TODO: Must be totally reimplemented (S.Zharko)
- static bool first = 1;
-
- ///Write Tracks and MC Tracks
-
- static int TrNumber = 0;
- std::fstream Tracks, McTracksCentr, McTracksIn, McTracksOut;
- if (first) {
- Tracks.open("tracks.dat", std::fstream::out);
- McTracksCentr.open("mctrackscentr.dat", std::fstream::out);
- McTracksIn.open("mctracksin.dat", std::fstream::out);
- McTracksOut.open("mctracksout.dat", std::fstream::out);
- first = 0;
- }
- else {
- Tracks.open("tracks.dat", std::fstream::out | std::fstream::app);
- McTracksCentr.open("mctrackscentr.dat", std::fstream::out | std::fstream::app);
- McTracksIn.open("mctracksin.dat", std::fstream::out | std::fstream::app);
- McTracksOut.open("mctracksout.dat", std::fstream::out | std::fstream::app);
- }
-
- for (auto RecTrack = fvRecoTracks.begin(); RecTrack != fvRecoTracks.end(); ++RecTrack) {
- if (RecTrack->IsGhost()) continue;
- const auto& mcTrk = fMCData.GetTrack(RecTrack->GetMatchedMCTrackIndex());
- if (!(mcTrk.IsPrimary())) continue;
-
- int NHits = (RecTrack->Hits).size();
- float x[20], y[20], z[20];
- int st[20];
- int jHit = 0;
- for (int iHit = 0; iHit < NHits; iHit++) {
- CbmL1HitDebugInfo& h = fvHitDebugInfo[RecTrack->Hits[iHit]];
- st[jHit] = h.iStation;
- if (h.ExtIndex < 0) {
- CbmMvdHit* MvdH = (CbmMvdHit*) fpMvdHits->At(-h.ExtIndex - 1);
- x[jHit] = MvdH->GetX();
- y[jHit] = MvdH->GetY();
- z[jHit] = MvdH->GetZ();
- jHit++;
- }
- else {
- CbmStsHit* StsH = (CbmStsHit*) fpStsHits->At(h.ExtIndex);
- x[jHit] = StsH->GetX();
- y[jHit] = StsH->GetY();
- z[jHit] = StsH->GetZ();
- jHit++;
- }
- }
-
- Tracks << TrNumber << " " << mcTrk.GetStartX() << " " << mcTrk.GetStartY() << " " << mcTrk.GetStartZ() << " "
- << mcTrk.GetPx() << " " << mcTrk.GetPy() << " " << mcTrk.GetPz() << " " << mcTrk.GetCharge() << " " << NHits
- << endl;
-
- float AngleXAxis = 0, AngleYAxis = 0;
- for (int i = 0; i < NHits; i++)
- Tracks << " " << st[i] << " " << x[i] << " " << y[i] << " " << z[i] << " " << AngleXAxis << " " << AngleYAxis
- << endl;
- Tracks << endl;
-
- int NMCPoints = mcTrk.GetNofPoints();
-
- McTracksIn << TrNumber << " " << mcTrk.GetStartX() << " " << mcTrk.GetStartY() << " " << mcTrk.GetStartZ() << " "
- << mcTrk.GetPx() << " " << mcTrk.GetPy() << " " << mcTrk.GetPz() << " " << mcTrk.GetCharge() << " "
- << NMCPoints << endl;
- McTracksOut << TrNumber << " " << mcTrk.GetStartX() << " " << mcTrk.GetStartY() << " " << mcTrk.GetStartZ() << " "
- << mcTrk.GetPx() << " " << mcTrk.GetPy() << " " << mcTrk.GetPz() << " " << mcTrk.GetCharge() << " "
- << NMCPoints << endl;
- McTracksCentr << TrNumber << " " << mcTrk.GetStartX() << " " << mcTrk.GetStartY() << " " << mcTrk.GetStartZ() << " "
- << mcTrk.GetPx() << " " << mcTrk.GetPy() << " " << mcTrk.GetPz() << " " << mcTrk.GetCharge() << " "
- << NMCPoints << endl;
-
- for (int iPoint : mcTrk.GetPointIndexes()) {
- const auto& mcPoint = fMCData.GetPoint(iPoint);
- McTracksIn << " " << mcPoint.GetStationId() << " " << mcPoint.GetXIn() << " " << mcPoint.GetYIn() << " "
- << mcPoint.GetZIn() << " " << mcPoint.GetPxIn() << " " << mcPoint.GetPyIn() << " " << mcPoint.GetPzIn()
- << endl;
- McTracksOut << " " << mcPoint.GetStationId() << " " << mcPoint.GetXOut() << " " << mcPoint.GetYOut() << " "
- << mcPoint.GetZOut() << " " << mcPoint.GetPxOut() << " " << mcPoint.GetPyOut() << " "
- << mcPoint.GetPzOut() << endl;
- McTracksCentr << " " << mcPoint.GetStationId() << " " << mcPoint.GetX() << " " << mcPoint.GetY() << " "
- << mcPoint.GetZ() << " " << mcPoint.GetPx() << " " << mcPoint.GetPy() << " " << mcPoint.GetPz()
- << endl;
- ;
- }
- McTracksIn << endl;
- McTracksOut << endl;
- McTracksCentr << endl;
-
- TrNumber++;
- }
-}
-
-
-// ---------------------------------------------------------------------------------------------------------------------
-//
-double CbmL1::boundedGaus(double sigma)
-{
- assert(sigma > 0. && std::isfinite(sigma));
- double x = 0.;
- do {
- x = gRandom->Gaus(0., sigma);
- } while (fabs(x) > 3.5 * sigma);
- return x;
-}
-
void CbmL1::DumpMaterialToFile(TString fileName)
{
TFile* f = new TFile(fileName, "RECREATE");
@@ -1318,9 +1043,7 @@ void CbmL1::DumpMaterialToFile(TString fileName)
}
// ---------------------------------------------------------------------------------------------------------------------
-// Target finder
//
-
void CbmL1::GetTargetInfo()
{
// Loop over all nodes till a node with name "target" is found
diff --git a/reco/L1/CbmL1.h b/reco/L1/CbmL1.h
index 077a6c24ce..a47d7442a6 100644
--- a/reco/L1/CbmL1.h
+++ b/reco/L1/CbmL1.h
@@ -36,12 +36,6 @@
#include "CbmL1MCPoint.h"
#include "CbmL1MCTrack.h"
#include "CbmL1Track.h"
-#include "CbmMCDataArray.h"
-#include "CbmMCEventList.h"
-#include "CbmMCTrack.h"
-#include "CbmMvdHit.h"
-#include "CbmMvdPoint.h"
-#include "CbmTimeSlice.h"
#include "FairDetector.h"
#include "FairRootManager.h"
#include "FairTask.h"
@@ -232,13 +226,6 @@ class CbmL1 : public FairTask {
/// \param filename Name of the output file name
void SetOutputMcTripletsTreeFilename(const char* filename) { fsMcTripletsOutputFilename = std::string(filename); }
- /// Sets flag: to correct input hits on MC or not
- /// \param flag: true - hits will be corrected on MC information
- void SetCorrectHitsOnMC(bool flag) { fIfCorrectHitsOnMC = flag; }
-
- /// Gets flag: to correct input hits on MC or not
- bool GetCorrectHitsOnMC() const { return fIfCorrectHitsOnMC; }
-
/// Gets vector of the extended QA hits
/// TODO: It is a temporary function, do not rely on it
const auto& GetQaHits() const { return fvHitDebugInfo; }
@@ -263,7 +250,7 @@ class CbmL1 : public FairTask {
void SetGlobalMode() { fTrackingMode = ca::Framework::TrackingMode::kGlobal; }
/// Sets misalignment of the detector
- void SetMisalignment(ca::EDetectorID det, float dx, float dy, float dt) { fvMisalignment[det] = {{dx, dy, dt}}; }
+ void SetMisalignment(ca::EDetectorID det, float dx, float dy, float dt) { fvMisalignment[det] = {dx, dy, dt}; }
void SetMisalignmentMvd(float dx, float dy, float dt) { SetMisalignment(ca::EDetectorID::kMvd, dx, dy, dt); }
void SetMisalignmentSts(float dx, float dy, float dt) { SetMisalignment(ca::EDetectorID::kSts, dx, dy, dt); }
@@ -353,7 +340,6 @@ class CbmL1 : public FairTask {
void FieldApproxCheck(); // Build histograms with difference between Field map and approximated field
void FieldIntegralCheck(); // Build 2D histogram: dependence of the field integral on phi and theta
- void InputPerformance(); // Build histograms about input data, like hit pulls, etc.
void TimeHist();
// ********************************
@@ -413,9 +399,6 @@ class CbmL1 : public FairTask {
void ReadSTAPPerfInputData();
- /// SIMD KF Banchmark service-functions
- void WriteSIMDKFData();
-
/// Gets a pointer to L1InitManager (for an access in run_reco.C)
ca::InitManager* GetInitManager() { return &fInitManager; }
@@ -425,7 +408,6 @@ class CbmL1 : public FairTask {
/// A helper for GetTargetInfo()
void FindTargetNode(TString& targetPath, TGeoNode*& targetNode);
-
void WriteHistosCurFile(TObject* obj);
static std::istream& eatwhite(std::istream& is); // skip spaces
@@ -448,10 +430,7 @@ class CbmL1 : public FairTask {
ca::InitManager fInitManager; ///< Tracking parameters data manager
std::shared_ptr<ca::DataManager> fpIODataManager = nullptr; ///< Input-output data manager
-
- //std::unique_ptr<CbmCaMCModule> fpMCModule = nullptr; ///< MC-module for tracking
-
- cbm::ca::DetIdArr_t<std::array<float, 3>> fvMisalignment{{{{{0.f, 0.f, 0.f}}}}}; ///< Misalignment
+ cbm::ca::DetIdArr_t<std::array<float, 3>> fvMisalignment{{0.}}; ///< Misalignment
public:
// ** Basic data members **
@@ -473,24 +452,10 @@ class CbmL1 : public FairTask {
double fTargetZ{1.e10}; ///< target position Z
int fNStations = 0; ///< number of total active detector stations
- int fNMvdStations = 0; ///< number of active MVD stations
- int fNStsStations = 0; ///< number of active STS stations
- int fNMuchStations = 0; ///< number of active MuCh stations
- int fNTrdStations = 0; ///< number of active TRD stations
- int fNTofStations = 0; ///< number of active TOF stations
-
- int fNStationsGeom = 0; ///< number of total detector stations
- int fNMvdStationsGeom = 0; ///< number of MVD stations
- int fNStsStationsGeom = 0; ///< number of STS stations;
- int fNMuchStationsGeom = 0; ///< number of MuCh stations;
- int fNTrdStationsGeom = 0; ///< number of TRD stations;
- int fNTofStationsGeom = 0; ///< number of TOF stations;
-
Int_t fPerformance = 0; ///< performance mode: 0 - w\o perf. 1 - L1-Efficiency definition. 2 - QA-Eff.definition
double fTrackingTime = 0.; ///< time of track finding procedure
-
/// Option to work with files for the standalone mode
/// - #0 standalone mode is not used
/// - #1 data for standalone mode is written to configuration file (currently does not work)
@@ -529,66 +494,6 @@ class CbmL1 : public FairTask {
bool fUseTRD = false; ///< if Trd data should be processed
bool fUseTOF = false; ///< if Tof data should be processed
- bool fIfCorrectHitsOnMC = false; ///< if correct input hits on MC information (debug)
-
- // ** Raw input data **
-
- CbmTimeSlice* fTimeSlice = nullptr; ///< Pointer to the TS object
-
- // Reconstructed hits input
- TClonesArray* fpMvdHits = nullptr; ///< Array of MVD hits ("MvdHit")
- TClonesArray* fpStsHits = nullptr; ///< Array of STS hits ("StsHit")
- TClonesArray* fpMuchPixelHits = nullptr; ///< Array of MuCh hits ("MuchPixelHit")
- TClonesArray* fpTrdHits = nullptr; ///< Array of TRD hits ("TrdHit")
- TClonesArray* fpTofHits = nullptr; ///< Array of TOF hits ("TofHit")
-
- // ** MC input **
-
- // General
- CbmMCEventList* fpMcEventList = nullptr; ///< MC event list (all)
- CbmMCDataArray* fpMcTracks = nullptr; ///< MC tracks list
- CbmMCDataObject* fpMcEventHeader = nullptr; ///< MC event header
-
- // MC-point arrays
- CbmMCDataArray* fpMvdPoints = nullptr; ///< MVD MC-points array
- CbmMCDataArray* fpStsPoints = nullptr; ///< STS MC-points array
- CbmMCDataArray* fpMuchPoints = nullptr; ///< MuCh MC-points array
- CbmMCDataArray* fpTrdPoints = nullptr; ///< TRD MC-points array
- CbmMCDataArray* fpTofPoints = nullptr; ///< TOF MC-points array
-
-
- // ** MC-reconstruction matching input **
-
- // Measured digis
- TClonesArray* fpMuchDigis = nullptr; ///< Array of MuCh digis
-
- // Reconstructed clusters
- TClonesArray* fpStsClusters = nullptr; ///< Array of STS clusters
- TClonesArray* fpMuchClusters = nullptr; ///< Array of MuCh clusters
-
- // Hit matches
- TClonesArray* fpStsHitMatches = nullptr; ///< Array of STS hit matches (NOTE: currently not used)
- TClonesArray* fpMvdHitMatches = nullptr; ///< Array of MVD hit matches
- TClonesArray* fpMuchHitMatches = nullptr; ///< Array of MuCh hit matches
- TClonesArray* fpTrdHitMatches = nullptr; ///< Array of TOF hit matches
- TClonesArray* fpTofHitMatches = nullptr; ///< Array of TRD hit matches
-
- // Digi and cluster matches
- TClonesArray* fpMvdDigiMatches = nullptr; ///< Array of MVD digi matches (NOTE: currently unused)
- TClonesArray* fpStsDigiMatches = nullptr; ///< Array of STS digi matches (NOTE: currently unused)
- TClonesArray* fpStsClusterMatches = nullptr; ///< Array of STS cluster matches
- TClonesArray* fpMuchDigiMatches = nullptr; ///< Array of MuCh digi matches (NOTE: currently unsused)
-
- ca::Vector<CbmL1MCPoint> fvMCPoints = {"CbmL1::fvMCPoints"}; ///< Container of MC points
-
-
- ca::Vector<CbmL1MCTrack> fvMCTracks = {"CbmL1::fvMCTracks"}; ///< Array of MC tracks
- ca::Vector<std::vector<int>> fvHitPointIndices = {
- "CbmL1::fvHitPointIndices"}; ///< Indices of MC points vs. hit index
- ca::Vector<int> fvHitBestPointIndices = {
- "CbmL1::fvHitBestPointIndices"}; ///< Indices of best MC points vs. hit index
-
-
int fEventNo = 0; ///< Current number of event/TS
int fNofRecoTracks = 0; ///< Total number of reconstructed tracks
@@ -603,10 +508,6 @@ class CbmL1 : public FairTask {
// indices of MCPoints in fvMCPoints, indexed by index of hit in algo->vHits array. According to StsMatch. Used for IdealResponce
// ca::Vector<int> vHitMCRef1;
// CbmMatch HitMatch;
- private:
- std::unordered_map<CbmL1LinkKey, int> fmMCPointsLinksMap = {}; /// Internal MC point index vs. link
- std::unordered_map<CbmL1LinkKey, int> fmMCTracksLinksMap = {}; /// Internal MC track index vs. link
-
cbm::ca::DetIdArr_t<std::set<int>> fvmDisabledStationIDs; /// Array of local indices of disabled tracking stations
diff --git a/reco/L1/CbmL1MCTrack.cxx b/reco/L1/CbmL1MCTrack.cxx
index 919f36c8c1..9aff8bdc4f 100644
--- a/reco/L1/CbmL1MCTrack.cxx
+++ b/reco/L1/CbmL1MCTrack.cxx
@@ -71,9 +71,8 @@ void CbmL1MCTrack::Init()
// get Hits
Hits.clear();
for (unsigned int iP = 0; iP < Points.size(); iP++) {
- CbmL1MCPoint* point = &(L1->fvMCPoints[Points[iP]]);
- for (unsigned int iH = 0; iH < point->hitIds.size(); iH++) {
- const int iih = point->hitIds[iH];
+ const auto& point = (L1->GetMCData().GetPoint(Points[iP]));
+ for (unsigned int iih : point.GetHitIndexes()) {
if (std::find(Hits.begin(), Hits.end(), iih) == Hits.end()) Hits.push_back_no_warning(iih);
}
}
@@ -92,8 +91,8 @@ void CbmL1MCTrack::CalculateMCCont()
nMCContStations = 0;
int istaold = -1, ncont = 0;
for (int ih = 0; ih < nPoints; ih++) {
- CbmL1MCPoint& h = L1->fvMCPoints[Points[ih]];
- int ista = h.iStation;
+ const auto& h = L1->GetMCData().GetPoint(Points[ih]);
+ int ista = h.GetStationId();
if (ista - istaold == 1)
ncont++;
else if (ista - istaold > 1) {
@@ -152,22 +151,22 @@ void CbmL1MCTrack::CalculateMaxNStaMC()
float lastz = -1;
int cur_maxNStaMC = 0, cur_maxNSensorMC = 0;
for (unsigned int iH = 0; iH < Points.size(); iH++) {
- CbmL1MCPoint& mcP = L1->fvMCPoints[Points[iH]];
- if (mcP.iStation == lastSta)
+ const auto& mcP = L1->GetMCData().GetPoint(Points[iH]);
+ if (mcP.GetStationId() == lastSta)
cur_maxNStaMC++;
else { // new station
if (cur_maxNStaMC > maxNStaMC) maxNStaMC = cur_maxNStaMC;
cur_maxNStaMC = 1;
- lastSta = mcP.iStation;
+ lastSta = mcP.GetStationId();
nMCStations++;
}
- if (mcP.z == lastz) // TODO: works incorrect because of different z
+ if (mcP.GetZ() == lastz) // TODO: works incorrect because of different z
cur_maxNSensorMC++;
else { // new z
if (cur_maxNSensorMC > maxNSensorMC) maxNSensorMC = cur_maxNSensorMC;
cur_maxNSensorMC = 1;
- lastz = mcP.z;
+ lastz = mcP.GetZ();
}
};
if (cur_maxNStaMC > maxNStaMC) maxNStaMC = cur_maxNStaMC;
@@ -198,7 +197,7 @@ void CbmL1MCTrack::CalculateIsReconstructable()
if (Points.size() > 0) {
isAdditional = f & (nHitContStations == nStations) & (nMCContStations == nStations) & (nMCStations == nStations)
- & (nHitContStations >= 3) & (L1->fvMCPoints[Points[0]].iStation == 0);
+ & (nHitContStations >= 3) & (L1->GetMCData().GetPoint(Points[0]).GetStationId() == 0);
isAdditional &= !isReconstructable;
}
else
diff --git a/reco/L1/CbmL1Performance.cxx b/reco/L1/CbmL1Performance.cxx
index 1b2d80943d..3de7372789 100644
--- a/reco/L1/CbmL1Performance.cxx
+++ b/reco/L1/CbmL1Performance.cxx
@@ -1748,491 +1748,6 @@ void CbmL1::FieldIntegralCheck()
gDirectory = curr;
} // void CbmL1::FieldApproxCheck()
-void CbmL1::InputPerformance()
-{
- // The input performance is currently evaluated by other CBM tasks
- // TODO: obsolete code, remove it
-
- // static TH1I *nStripFHits, *nStripBHits, *nStripFMC, *nStripBMC;
-
- static TH1F *resXsts, *resYsts, *resTsts, *resXmvd, *resYmvd /*, *pullZ*/;
- static TH1F *pullXsts, *pullYsts, *pullTsts, *pullXmvd, *pullYmvd /*, *pullZ*/;
- static TH1F *pullXmuch, *pullYmuch, *pullTmuch, *resXmuch, *resYmuch, *resTmuch /*, *pullZ*/;
- static TH1F *pullXtrd, *pullYtrd, *pullTtrd, *resXtrd, *resYtrd, *resTtrd /*, *pullZ*/;
- static TH1F *pullXtof, *pullYtof, *pullTtof, *resXtof, *resYtof, *resTtof /*, *pullZ*/;
- static TH1F* trdMCPointsZ = nullptr;
-
- static bool first_call = 1;
- if (first_call) {
- first_call = 0;
-
- TDirectory* currentDir = gDirectory;
- gDirectory = fHistoDir;
- gDirectory->cd("Input");
- gDirectory->mkdir("STS");
- gDirectory->cd("STS");
-
- // nStripFHits = new TH1I("nHits_f", "NHits On Front Strip", 10, 0, 10);
- // nStripBHits = new TH1I("nHits_b", "NHits On Back Strip", 10, 0, 10);
- // nStripFMC = new TH1I("nMC_f", "N MC Points On Front Strip", 10, 0, 10);
- // nStripBMC = new TH1I("nMC_b", "N MC Points On Back Strip", 10, 0, 10);
-
- pullXsts = new TH1F("Px_sts", "STS Pull x", 100, -5, 5);
- pullYsts = new TH1F("Py_sts", "STS Pull y", 100, -5, 5);
- pullTsts = new TH1F("Pt_sts", "STS Pull t", 100, -5, 5);
-
- resXsts = new TH1F("x_sts", "STS Residual x", 100, -50, 50);
- resYsts = new TH1F("y_sts", "STS Residual y", 100, -500, 500);
- resTsts = new TH1F("t_sts", "STS Residual t", 100, -20, 20);
-
- gDirectory->cd("..");
- gDirectory->mkdir("MVD");
- gDirectory->cd("MVD");
-
- pullXmvd = new TH1F("Px_mvd", "MVD Pull x", 100, -5, 5);
- pullYmvd = new TH1F("Py_mvd", "MVD Pull y", 100, -5, 5);
-
- resXmvd = new TH1F("x_mvd", "MVD Residual x", 100, -50, 50);
- resYmvd = new TH1F("y_mvd", "MVD Residual y", 100, -50, 50);
-
- TH1* histo;
- histo = resXsts;
- histo->GetXaxis()->SetTitle("Residual x, um");
- histo = resYsts;
- histo->GetXaxis()->SetTitle("Residual y, um");
- histo = resTsts;
- histo->GetXaxis()->SetTitle("Residual t, ns");
- histo = resXmvd;
- histo->GetXaxis()->SetTitle("Residual x, um");
- histo = resYmvd;
- histo->GetXaxis()->SetTitle("Residual y, um");
- histo = pullXsts;
- histo->GetXaxis()->SetTitle("Pull x");
- histo = pullYsts;
- histo->GetXaxis()->SetTitle("Pull y");
- histo = pullTsts;
- histo->GetXaxis()->SetTitle("Pull t");
- histo = pullXmvd;
- histo->GetXaxis()->SetTitle("Pull x");
- histo = pullYmvd;
- histo->GetXaxis()->SetTitle("Pull y");
-
-
- gDirectory->cd("..");
- gDirectory->mkdir("MUCH");
- gDirectory->cd("MUCH");
-
- pullXmuch = new TH1F("Px_much", "MUCH Pull x", 100, -10, 10);
- pullYmuch = new TH1F("Py_much", "MUCH Pull y", 100, -10, 10);
- pullTmuch = new TH1F("Pt_much", "MUCH Pull t", 100, -10, 10);
-
- resXmuch = new TH1F("x_much", "MUCH Residual x", 100, -100000, 100000);
- resYmuch = new TH1F("y_much", "MUCH Residual y", 100, -100000, 100000);
- resTmuch = new TH1F("t_much", "MUCH Residual t", 100, -40, 40);
-
-
- histo = resXmuch;
- histo->GetXaxis()->SetTitle("Residual x, um");
- histo = resYmuch;
- histo->GetXaxis()->SetTitle("Residual y, um");
- histo = resTmuch;
- histo->GetXaxis()->SetTitle("Residual t, ns");
- histo = pullXmuch;
- histo->GetXaxis()->SetTitle("Pull x");
- histo = pullYmuch;
- histo->GetXaxis()->SetTitle("Pull y");
- histo = pullTmuch;
- histo->GetXaxis()->SetTitle("Pull t");
-
- gDirectory->cd("..");
- gDirectory->mkdir("TRD");
- gDirectory->cd("TRD");
-
- trdMCPointsZ = new TH1F("trdMCPointsZ", "trdMCPointsZ", 1000, 400., 700.);
-
- pullXtrd = new TH1F("Px_trd", "TRD Pull x", 100, -5, 5);
- pullYtrd = new TH1F("Py_trd", "TRD Pull y", 100, -5, 5);
- pullTtrd = new TH1F("Pt_trd", "TRD Pull t", 100, -5, 5);
-
- resXtrd = new TH1F("x_trd", "TRD Residual x", 100, -10000, 10000);
- resYtrd = new TH1F("y_trd", "TRD Residual y", 100, -10000, 10000);
- resTtrd = new TH1F("t_trd", "TRD Residual t", 100, -40, 40);
-
-
- histo = resXtrd;
- histo->GetXaxis()->SetTitle("Residual x, um");
- histo = resYtrd;
- histo->GetXaxis()->SetTitle("Residual y, um");
- histo = resTtrd;
- histo->GetXaxis()->SetTitle("Residual t, ns");
- histo = pullXtrd;
- histo->GetXaxis()->SetTitle("Pull x");
- histo = pullYtrd;
- histo->GetXaxis()->SetTitle("Pull y");
- histo = pullTtrd;
- histo->GetXaxis()->SetTitle("Pull t");
-
- gDirectory->cd("..");
- gDirectory->mkdir("TOF");
- gDirectory->cd("TOF");
-
- pullXtof = new TH1F("Px_tof", "TOF Pull x", 100, -5, 5);
- pullYtof = new TH1F("Py_tof", "TOF Pull y", 100, -5, 5);
- pullTtof = new TH1F("Pt_tof", "TOF Pull t", 100, -5, 5);
-
- resXtof = new TH1F("x_tof", "TOF Residual x", 100, -50000, 50000);
- resYtof = new TH1F("y_tof", "TOF Residual y", 100, -50000, 50000);
- resTtof = new TH1F("t_tof", "TOF Residual t", 100, -0.4, 0.4);
-
-
- histo = resXtof;
- histo->GetXaxis()->SetTitle("Residual x, um");
- histo = resYtof;
- histo->GetXaxis()->SetTitle("Residual y, um");
- histo = resTtof;
- histo->GetXaxis()->SetTitle("Residual t, ns");
- histo = pullXtof;
- histo->GetXaxis()->SetTitle("Pull x");
- histo = pullYtof;
- histo->GetXaxis()->SetTitle("Pull y");
- histo = pullTtof;
- histo->GetXaxis()->SetTitle("Pull t");
-
- gDirectory = currentDir;
- } // first_call
-
- // std::map<unsigned int, unsigned int> stripFToNHitMap,stripBToNHitMap;
- // std::map<unsigned int, unsigned int> stripFToNMCMap,stripBToNMCMap;
-
- map<unsigned int, unsigned int>::iterator it;
-
- if (fpStsHits && fpStsHitMatches && fpStsPoints) {
- for (int iH = 0; iH < fpStsHits->GetEntriesFast(); iH++) {
-
- const CbmStsHit* sh = dynamic_cast<CbmStsHit*>(fpStsHits->At(iH));
-
- // int iMCPoint = -1;
- CbmLink link;
- CbmStsPoint* pt = nullptr;
-
- if (fpStsClusterMatches) {
- const CbmMatch* frontClusterMatch =
- static_cast<const CbmMatch*>(fpStsClusterMatches->At(sh->GetFrontClusterId()));
- const CbmMatch* backClusterMatch =
- static_cast<const CbmMatch*>(fpStsClusterMatches->At(sh->GetBackClusterId()));
- CbmMatch stsHitMatch;
-
- for (Int_t iFrontLink = 0; iFrontLink < frontClusterMatch->GetNofLinks(); iFrontLink++) {
- const CbmLink& frontLink = frontClusterMatch->GetLink(iFrontLink);
-
- for (Int_t iBackLink = 0; iBackLink < backClusterMatch->GetNofLinks(); iBackLink++) {
- const CbmLink& backLink = backClusterMatch->GetLink(iBackLink);
- if (frontLink == backLink) {
- stsHitMatch.AddLink(frontLink);
- stsHitMatch.AddLink(backLink);
- }
- }
- }
-
- if (stsHitMatch.GetNofLinks() > 0) {
- Float_t bestWeight = 0.f;
- for (Int_t iLink = 0; iLink < stsHitMatch.GetNofLinks(); iLink++) {
- if (stsHitMatch.GetLink(iLink).GetWeight() > bestWeight) {
- link = stsHitMatch.GetLink(iLink);
- if (link.GetIndex() < 0) break;
- bestWeight = link.GetWeight();
- Int_t iFile = link.GetFile();
- Int_t iEvent = link.GetEntry();
- pt = (CbmStsPoint*) fpStsPoints->Get(iFile, iEvent, link.GetIndex());
- }
- }
- }
-
- if (pt == nullptr) continue;
-
- double mcTime = pt->GetTime();
-
- mcTime += fpMcEventList->GetEventTime(link.GetEntry(), link.GetFile());
-
- // hit pulls and residuals
-
- TVector3 hitPos, mcPos, hitErr;
- sh->Position(hitPos);
- sh->PositionError(hitErr);
-
- double t = (hitPos.Z() - pt->GetZIn()) / (pt->GetZOut() - pt->GetZIn());
- mcPos.SetX(pt->GetXIn() + t * (pt->GetXOut() - pt->GetXIn()));
- mcPos.SetY(pt->GetYIn() + t * (pt->GetYOut() - pt->GetYIn()));
- mcPos.SetZ(hitPos.Z());
-
- { // qa errors
- if (sh->GetDx() > 1.e-8) pullXsts->Fill((hitPos.X() - mcPos.X()) / sh->GetDx());
- if (sh->GetDy() > 1.e-8) pullYsts->Fill((hitPos.Y() - mcPos.Y()) / sh->GetDy());
- pullTsts->Fill((sh->GetTime() - mcTime) / sh->GetTimeError());
- }
-
- resXsts->Fill((hitPos.X() - mcPos.X()) * 10 * 1000);
- resYsts->Fill((hitPos.Y() - mcPos.Y()) * 10 * 1000);
- resTsts->Fill((sh->GetTime() - mcTime));
- }
- }
- } // sts
-
-
- if (fpMvdHits && fpMvdHitMatches && fpMvdPoints) {
- Int_t nEnt = fpMvdHits->GetEntriesFast();
- for (int j = 0; j < nEnt; j++) {
-
- CbmMvdHit* sh = dynamic_cast<CbmMvdHit*>(fpMvdHits->At(j));
- CbmMatch* hm = dynamic_cast<CbmMatch*>(fpMvdHitMatches->At(j));
-
- CbmMvdPoint* pt = nullptr;
- {
- float mcWeight = -1.f;
- for (int iLink = 0; iLink < hm->GetNofLinks(); iLink++) {
- const CbmLink& link = hm->GetLink(iLink);
- if (link.GetIndex() < 0) break;
- if (link.GetWeight() < mcWeight) continue;
- mcWeight = link.GetWeight();
- pt = dynamic_cast<CbmMvdPoint*>(fpMvdPoints->Get(&link));
- }
- }
- if (!pt) continue;
-
- // hit pulls and residuals
-
- TVector3 hitPos, mcPos, hitErr;
- sh->Position(hitPos);
- sh->PositionError(hitErr);
-
- mcPos.SetX((pt->GetX() + pt->GetXOut()) / 2.);
- mcPos.SetY((pt->GetY() + pt->GetYOut()) / 2.);
- mcPos.SetZ(hitPos.Z());
-
- if (hitErr.X() != 0) pullXmvd->Fill((hitPos.X() - mcPos.X()) / hitErr.X());
- if (hitErr.Y() != 0) pullYmvd->Fill((hitPos.Y() - mcPos.Y()) / hitErr.Y());
-
- resXmvd->Fill((hitPos.X() - mcPos.X()) * 10 * 1000);
- resYmvd->Fill((hitPos.Y() - mcPos.Y()) * 10 * 1000);
- }
- } // mvd
-
-
- if (fpMuchPixelHits && fpMuchHitMatches && fpMuchPoints) {
- for (int iH = 0; iH < fpMuchPixelHits->GetEntriesFast(); iH++) {
-
- const CbmMuchPixelHit* sh = dynamic_cast<CbmMuchPixelHit*>(fpMuchPixelHits->At(iH));
- const CbmMatch* hm = dynamic_cast<CbmMatch*>(fpMuchHitMatches->At(iH));
-
- if (hm->GetNofLinks() == 0) continue;
- Float_t bestWeight = 0.f;
- Float_t totalWeight = 0.f;
- int iMCPoint = -1;
- CbmLink link;
-
- for (int iLink = 0; iLink < hm->GetNofLinks(); iLink++) {
- totalWeight += hm->GetLink(iLink).GetWeight();
- if (hm->GetLink(iLink).GetWeight() > bestWeight) {
- if (hm->GetLink(iLink).GetIndex() < 0) break;
- bestWeight = hm->GetLink(iLink).GetWeight();
- iMCPoint = hm->GetLink(iLink).GetIndex();
- link = hm->GetLink(iLink);
- }
- }
- if (bestWeight / totalWeight < 0.7 || iMCPoint < 0) continue;
-
- CbmMuchPoint* pt = (CbmMuchPoint*) fpMuchPoints->Get(link.GetFile(), link.GetEntry(), link.GetIndex());
- double mcTime = pt->GetTime();
- mcTime += fpMcEventList->GetEventTime(link.GetEntry(), link.GetFile());
- // mcTime+=20;
-
- // hit pulls and residuals
-
-
- TVector3 hitPos, mcPos, hitErr;
- sh->Position(hitPos);
- sh->PositionError(hitErr);
-
- // pt->Position(mcPos); // this is wrong!
- // mcPos.SetX( pt->GetX( hitPos.Z() ) );
- // mcPos.SetY( pt->GetY( hitPos.Z() ) );
- // mcPos.SetZ( hitPos.Z() );
-
- mcPos.SetX(0.5 * (pt->GetXIn() + pt->GetXOut()));
- mcPos.SetY(0.5 * (pt->GetYIn() + pt->GetYOut()));
- mcPos.SetZ(hitPos.Z());
-
- { // qa errors
- if (sh->GetDx() > 1.e-8) pullXmuch->Fill((sh->GetX() - mcPos.X()) / sh->GetDx());
- if (sh->GetDy() > 1.e-8) pullYmuch->Fill((sh->GetY() - mcPos.Y()) / sh->GetDy());
- pullTmuch->Fill((sh->GetTime() - mcTime) / sh->GetTimeError());
- }
-
- resXmuch->Fill((sh->GetX() - mcPos.X()) * 10 * 1000);
- resYmuch->Fill((sh->GetY() - mcPos.Y()) * 10 * 1000);
- resTmuch->Fill((sh->GetTime() - mcTime));
- }
- } // much
-
-
- if (fpTrdHits && fpTrdHitMatches && fpTrdPoints) {
- for (int iH = 0; iH < fpTrdHits->GetEntriesFast(); iH++) {
-
- const CbmTrdHit* sh = dynamic_cast<CbmTrdHit*>(fpTrdHits->At(iH));
- const CbmMatch* hm = dynamic_cast<CbmMatch*>(fpTrdHitMatches->At(iH));
-
- if (hm->GetNofLinks() == 0) continue;
- if (hm->GetNofLinks() != 1) continue; // only check single-linked hits
-
- Float_t bestWeight = 0.f;
- Float_t totalWeight = 0.f;
- int iMCPoint = -1;
- CbmLink link;
-
- for (int iLink = 0; iLink < hm->GetNofLinks(); iLink++) {
- totalWeight += hm->GetLink(iLink).GetWeight();
- if (hm->GetLink(iLink).GetWeight() > bestWeight) {
- if (hm->GetLink(iLink).GetIndex() < 0) {
- iMCPoint = -1;
- break;
- }
- bestWeight = hm->GetLink(iLink).GetWeight();
- iMCPoint = hm->GetLink(iLink).GetIndex();
- link = hm->GetLink(iLink);
- }
- }
- if (bestWeight / totalWeight < 0.7 || iMCPoint < 0) continue;
-
- CbmTrdPoint* pt = (CbmTrdPoint*) fpTrdPoints->Get(link.GetFile(), link.GetEntry(), link.GetIndex());
- double mcTime = pt->GetTime();
-
- mcTime += fpMcEventList->GetEventTime(link.GetEntry(), link.GetFile());
-
- // hit pulls and residuals
- // if ((sh->GetPlaneId()) == 0) continue;
- // if ((sh->GetPlaneId()) == 2) continue;
- // if ((sh->GetPlaneId()) == 4) continue;
-
- TVector3 hitPos, mcPos, hitErr;
- sh->Position(hitPos);
- sh->PositionError(hitErr);
-
- // pt->Position(mcPos); // this is wrong!
- mcPos.SetX((pt->GetXIn() + pt->GetXOut()) / 2.);
- mcPos.SetY((pt->GetYIn() + pt->GetYOut()) / 2.);
- mcPos.SetZ(hitPos.Z());
-
- { // qa errors
- if (sh->GetDx() > 1.e-8) pullXtrd->Fill((sh->GetX() - mcPos.X()) / sh->GetDx());
- if (sh->GetDy() > 1.e-8) pullYtrd->Fill((sh->GetY() - mcPos.Y()) / sh->GetDy());
- pullTtrd->Fill((sh->GetTime() - mcTime) / sh->GetTimeError());
- }
-
- resXtrd->Fill((sh->GetX() - mcPos.X()) * 10 * 1000);
- resYtrd->Fill((sh->GetY() - mcPos.Y()) * 10 * 1000);
- resTtrd->Fill((sh->GetTime() - mcTime));
- trdMCPointsZ->Fill(mcPos.Z());
- }
- } // much
-
-
- if (fpTofHits && fpTofHitMatches && fpTofPoints) {
- for (int iH = 0; iH < fpTofHits->GetEntriesFast(); iH++) {
-
- const CbmTofHit* sh = dynamic_cast<CbmTofHit*>(fpTofHits->At(iH));
- const CbmMatch* hm = dynamic_cast<CbmMatch*>(fpTofHitMatches->At(iH));
-
- if (hm->GetNofLinks() == 0) continue;
- Float_t bestWeight = 0.f;
- //Float_t totalWeight = 0.f;
- int iMCPoint = -1;
- CbmLink link;
-
- for (int iLink = 0; iLink < hm->GetNofLinks(); iLink++) {
- //totalWeight += hm->GetLink(iLink).GetWeight();
- if (hm->GetLink(iLink).GetWeight() > bestWeight) {
- if (hm->GetLink(iLink).GetIndex() < 0) {
- iMCPoint = -1;
- break;
- };
- bestWeight = hm->GetLink(iLink).GetWeight();
- iMCPoint = hm->GetLink(iLink).GetIndex();
- link = hm->GetLink(iLink);
- }
- }
-
- // TODO: Unify the cut and the whole code block for different detectors (S.Zharko)
- // if (bestWeight / totalWeight < 0.7 || iMCPoint < 0) continue;
- if (iMCPoint < 0) continue;
-
- CbmTofPoint* pt = (CbmTofPoint*) fpTofPoints->Get(link.GetFile(), link.GetEntry(), link.GetIndex());
- double mcTime = pt->GetTime();
-
- mcTime += fpMcEventList->GetEventTime(link.GetEntry(), link.GetFile());
-
- // hit pulls and residuals
-
-
- TVector3 hitPos, mcPos, hitErr;
- sh->Position(hitPos);
- sh->PositionError(hitErr);
-
- // pt->Position(mcPos); // this is wrong!
- mcPos.SetX((pt->GetX()));
- mcPos.SetY((pt->GetY()));
- mcPos.SetZ(hitPos.Z());
-
- { // qa errors
- if (sh->GetDx() > 1.e-8) pullXtof->Fill((sh->GetX() - mcPos.X()) / sh->GetDx());
- if (sh->GetDy() > 1.e-8) pullYtof->Fill((sh->GetY() - mcPos.Y()) / sh->GetDy());
- pullTtof->Fill((sh->GetTime() - mcTime) / sh->GetTimeError());
- }
-
- resXtof->Fill((sh->GetX() - mcPos.X()) * 10 * 1000);
- resYtof->Fill((sh->GetY() - mcPos.Y()) * 10 * 1000);
- resTtof->Fill((sh->GetTime() - mcTime));
- }
- } // much
-
- // for (it = stripFToNHitMap.begin(); it != stripFToNHitMap.end(); it++){
- // nStripFHits->Fill(it->second);
- // }
- // for (it = stripBToNHitMap.begin(); it != stripBToNHitMap.end(); it++){
- // nStripBHits->Fill(it->second);
- // }
- // for (it = stripFToNMCMap.begin(); it != stripFToNMCMap.end(); it++){
- // nStripFMC->Fill(it->second);
- // }
- // for (it = stripBToNMCMap.begin(); it != stripBToNMCMap.end(); it++){
- // nStripBMC->Fill(it->second);
- // }
-
- // strips Not ended
- // if( listCbmStsDigi ){
- // Int_t nEnt = listCbmStsDigi->GetEntriesFast();
- // for (int j=0; j < nEnt; j++ ){
- // CbmStsDigi *digi = (CbmStsDigi*) listCbmStsDigi->At(j);
- // // = sh->GetNLinks(0);
- // // find position of mc point
- // FairLink link = digi->GetLink(0);
- // int iMCPoint = link.GetIndex();
- // CbmStsPoint *mcPoint = (CbmStsPoint*) listStsPts->At(iMCPoint);
- // TVector3 mcPos;
- // if (digi->GetSide() == 0)
- // mcPoint->PositionIn(mcPos);
- // else
- // mcPoint->PositionOut(mcPos);
- //
- // CbmStsStation_old *sta = StsDigi.GetStation(digi->GetStationNr());
- // CbmStsSector* sector = sta->GetSector(digi->GetSectorNr());
- // digi->GetChannelNr();
- //
- // }
- // } // listCbmStsDigi
-
-
-} // void CbmL1::InputPerformance()
-
// ---------------------------------------------------------------------------------------------------------------------
//
void CbmL1::DumpMCTripletsToTree()
diff --git a/reco/L1/L1Algo/utils/L1AlgoPulls.cxx b/reco/L1/L1Algo/utils/L1AlgoPulls.cxx
index 9c4064ce91..2eb4ca5ce5 100644
--- a/reco/L1/L1Algo/utils/L1AlgoPulls.cxx
+++ b/reco/L1/L1Algo/utils/L1AlgoPulls.cxx
@@ -132,7 +132,7 @@ void L1AlgoPulls::AddOne(TrackParamV& T_, int i, ca::HitIndex_t ih)
// mc data
int iMCP = fL1->GetQaHits()[ih].GetBestMcPointId();
if (iMCP < 0) return;
- CbmL1MCPoint& mcP = fL1->fvMCPoints[iMCP];
+ const auto& mcP = fL1->GetMCData().GetPoint(iMCP);
TL1TrackParameters mc(mcP);
// fill residuals
diff --git a/reco/L1/L1Algo/utils/L1AlgoPulls.h b/reco/L1/L1Algo/utils/L1AlgoPulls.h
index c6fac2cc33..1014b9c775 100644
--- a/reco/L1/L1Algo/utils/L1AlgoPulls.h
+++ b/reco/L1/L1Algo/utils/L1AlgoPulls.h
@@ -43,6 +43,13 @@ struct TL1TrackParameters {
TL1TrackParameters(CbmL1MCPoint& T) : x(T.x), y(T.y), tx(T.px / T.pz), ty(T.py / T.pz), qp(T.q / T.p){};
+ TL1TrackParameters(const cbm::ca::tools::MCPoint& T)
+ : x(T.GetX())
+ , y(T.GetY())
+ , tx(T.GetTx())
+ , ty(T.GetTy())
+ , qp(T.GetQp()){};
+
double operator[](int i)
{
switch (i) {
diff --git a/reco/L1/qa/CbmCaOutputQa.cxx b/reco/L1/qa/CbmCaOutputQa.cxx
index 0db50630c9..ea25ca2d1e 100644
--- a/reco/L1/qa/CbmCaOutputQa.cxx
+++ b/reco/L1/qa/CbmCaOutputQa.cxx
@@ -938,8 +938,6 @@ InitStatus OutputQa::InitTimeSlice()
nHits = fvHits.size();
nRecoTracks = fvRecoTracks.size();
- LOG(info) << "Reco hits: " << nHits << ", reco tracks: " << nRecoTracks;
-
fMonitor.IncrementCounter(EMonitorKey::kTrack, nRecoTracks);
fMonitor.IncrementCounter(EMonitorKey::kHit, nHits);
@@ -958,58 +956,6 @@ InitStatus OutputQa::InitTimeSlice()
DrawEvent();
}
}
-
- // ****** DEBUG: BEGIN
- if constexpr (0) {
- using std::setw;
- {
- LOG(info) << "---------- Reco track sample";
- int id = 0;
- if (fvRecoTracks.size()) {
- LOG(info) << setw(4) << "No." << ' ' << fvRecoTracks[0].ToString(3, true);
- }
- for (const auto& trk : fvRecoTracks) {
- if (trk.IsGhost()) {
- LOG(info) << "I AM A GHOST! BOO!";
- }
- LOG(info) << setw(4) << (id++) << ' ' << trk.ToString(3);
- }
- }
- {
- LOG(info) << "---------- Reco hit sample";
- int id = 0;
- if (fvHits.size()) {
- LOG(info) << setw(4) << "No." << ' ' << fvHits[0].ToString(3, true);
- }
- for (const auto& h : fvHits) {
- LOG(info) << setw(4) << (id++) << ' ' << h.ToString(3);
- }
- }
- if (IsMCUsed()) {
- LOG(info) << "---------- MC tracks sample";
- int id = 0;
- if (fMCData.GetTrackContainer().size()) {
- LOG(info) << setw(4) << "No." << ' ' << fMCData.GetTrack(0).ToString(3, true);
- }
- for (const auto& trk : fMCData.GetTrackContainer()) {
- if (trk.GetNofHits() || trk.GetNofPoints()) {
- LOG(info) << setw(4) << (id++) << ' ' << trk.ToString(3);
- }
- }
- }
- if (IsMCUsed()) {
- LOG(info) << "---------- MC points sample";
- int id = 0;
- if (fMCData.GetPointContainer().size()) {
- LOG(info) << setw(4) << "No." << ' ' << fMCData.GetPoint(0).ToString(3, true);
- }
- for (const auto& p : fMCData.GetPointContainer()) {
- LOG(info) << setw(4) << (id++) << ' ' << p.ToString(3);
- }
- }
- }
- // ****** DEBUG: END
-
LOG_IF(info, fVerbose > 1) << fName << ": Data sample consists of " << nHits << " hits, " << nRecoTracks
<< " reco tracks, " << nMCTracks << " MC tracks, " << nMCPoints << " MC points";
--
GitLab