diff --git a/macro/run/run_qa.C b/macro/run/run_qa.C
index 1754454a3cd5ec45a9c03ee91f3646c42feb9fa0..be7fb5b0174c7e5539e5fb42fdc576b97870e0e2 100644
--- a/macro/run/run_qa.C
+++ b/macro/run/run_qa.C
@@ -43,6 +43,8 @@ void run_qa(TString dataTra = "data/sis100_muon_jpsi_test", TString dataRaw = "d
TString dataTra2 = "", TString dataTra3 = "")
{
+ gROOT->SetBatch(kTRUE);
+
// ========================================================================
// Adjust this part according to your requirements
diff --git a/reco/KF/CbmKF.cxx b/reco/KF/CbmKF.cxx
index 8ce12f2c8c463c556bd013c5b90677d92927eb0c..8d8c2f8c31f72d85492104c3edca0634746ba13a 100644
--- a/reco/KF/CbmKF.cxx
+++ b/reco/KF/CbmKF.cxx
@@ -262,65 +262,78 @@ void CbmKF::GetTargetInfo()
// 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
- TGeoNode* topNode = gGeoManager->GetTopNode();
-
- TObjArray* nodes = topNode->GetNodes();
-
- loop_over_nodes(nodes);
CbmKFTube target {};
target.ID = -111;
target.F = 1.;
- if (fTarget) {
- TGeoMatrix* matrix = fTarget->GetMatrix();
- const Double_t* translation = matrix->GetTranslation();
- target.x = translation[0];
- target.y = translation[1];
- target.z = translation[2];
+ TString targetPath;
+ TGeoNode* targetNode {nullptr};
+ FindTargetNode(targetPath, targetNode);
- TGeoVolume* volume = fTarget->GetVolume();
+ if (!targetNode) { LOG(fatal) << "Could not find the target."; }
- TGeoShape* shape = volume->GetShape();
- if (shape->TestShapeBit(TGeoShape::kGeoTube)) {
- target.r = static_cast<TGeoTube*>(shape)->GetRmin();
- target.R = static_cast<TGeoTube*>(shape)->GetRmax();
- target.dz = 2. * static_cast<TGeoTube*>(shape)->GetDz();
- }
- else {
- LOG(fatal) << "Only a target of a tube shape is supported";
- }
+ 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);
+ target.x = global[0];
+ target.y = global[1];
+ target.z = global[2];
- TGeoMaterial* material = volume->GetMaterial();
- Double_t radlength = material->GetRadLen();
- target.RadLength = radlength;
- target.Fe = 0.02145;
+ if (fVerbose) {
+ cout << "KALMAN FILTER : === READ TARGET MATERIAL ===" << endl;
+ cout << " found targed \"" << targetPath << "\" at ( " << target.x << " " << target.y << " " << target.z << " ) "
+ << endl;
+ }
- target.rr = target.r * target.r;
- target.RR = target.R * target.R;
- target.ZThickness = target.dz;
- target.ZReference = target.z;
+ TGeoVolume* volume = targetNode->GetVolume();
- vTargets.push_back(target);
- LOG(info) << "Target info: " << target.KFInfo();
+ TGeoShape* shape = volume->GetShape();
+ if (shape->TestShapeBit(TGeoShape::kGeoTube)) {
+ target.r = static_cast<TGeoTube*>(shape)->GetRmin();
+ target.R = static_cast<TGeoTube*>(shape)->GetRmax();
+ target.dz = 2. * static_cast<TGeoTube*>(shape)->GetDz();
}
else {
- LOG(fatal) << "Could not find the target.";
+ LOG(fatal) << "Only a target of a tube shape is supported";
}
+
+ TGeoMaterial* material = volume->GetMaterial();
+ Double_t radlength = material->GetRadLen();
+ target.RadLength = radlength;
+ target.Fe = 0.02145;
+
+ target.rr = target.r * target.r;
+ target.RR = target.R * target.R;
+ target.ZThickness = target.dz;
+ target.ZReference = target.z;
+
+ vTargets.push_back(target);
+ LOG(info) << "Target info: " << target.KFInfo();
}
-void CbmKF::loop_over_nodes(TObjArray* nodes)
+void CbmKF::FindTargetNode(TString& targetPath, TGeoNode*& targetNode)
{
- for (Int_t iNode = 0; iNode < nodes->GetEntriesFast(); iNode++) {
- TGeoNode* node = static_cast<TGeoNode*>(nodes->At(iNode));
- TString nodeName = node->GetName();
- if (nodeName.Contains("target")) {
- fTarget = node;
- break;
+ 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;
}
- TObjArray* subnodes = node->GetNodes();
- if (nullptr != subnodes) { loop_over_nodes(subnodes); }
}
+ targetPath = "";
+ targetNode = nullptr;
}
Int_t CbmKF::GetMaterialIndex(Int_t uid)
diff --git a/reco/KF/CbmKF.h b/reco/KF/CbmKF.h
index 7b98f9825f44125a4e2142109bba511c6b7a0ef6..c28d7b9981a8123e2bf1792c20e5376409871a33 100644
--- a/reco/KF/CbmKF.h
+++ b/reco/KF/CbmKF.h
@@ -107,9 +107,8 @@ private:
Int_t ReadTube(CbmKFTube& tube, FairGeoNode* node);
CbmKFMaterial* ReadPassive(FairGeoNode* node);
void GetTargetInfo();
- void loop_over_nodes(TObjArray* nodes);
- TGeoNode* fTarget {nullptr};
+ void FindTargetNode(TString& targetPath, TGeoNode*& targetNode);
private:
CbmKF(const CbmKF&);
diff --git a/reco/L1/CbmL1.cxx b/reco/L1/CbmL1.cxx
index 58a45a359b2639e36134f8c747fa740093cfade0..39d49953b89573395f24f3d703d0d73fc7a76479 100644
--- a/reco/L1/CbmL1.cxx
+++ b/reco/L1/CbmL1.cxx
@@ -31,6 +31,7 @@
#include "CbmMvdStationPar.h"
// TODO: include of CbmSetup.h creates problems on Mac
// #include "CbmSetup.h"
+#include "CbmMCDataObject.h"
#include "CbmStsFindTracks.h"
#include "CbmStsHit.h"
#include "CbmStsParSetModule.h"
@@ -285,9 +286,13 @@ InitStatus CbmL1::Init()
fStsPoints = mcManager->InitBranch("StsPoint");
+ fMcEventHeader = mcManager->GetObject("MCEventHeader.");
+
fMCTracks = mcManager->InitBranch("MCTrack");
+
if (NULL == fStsPoints) LOG(fatal) << GetName() << ": No StsPoint data!";
if (NULL == fMCTracks) LOG(fatal) << GetName() << ": No MCTrack data!";
+ if (NULL == fMcEventHeader) LOG(fatal) << GetName() << ": No MC event header data!";
if (!fLegacyEventMode) {
fEventList = (CbmMCEventList*) fManger->GetObject("MCEventList.");
@@ -353,14 +358,22 @@ InitStatus CbmL1::Init()
L1Vector<fscal> geo("geo");
geo.reserve(10000);
- for (int i = 0; i < 3; i++) {
- Double_t point[3] = {0, 0, 2.5 * i};
- Double_t B[3] = {0, 0, 0};
- if (CbmKF::Instance()->GetMagneticField()) CbmKF::Instance()->GetMagneticField()->GetFieldValue(point, B);
- geo.push_back(2.5 * i);
- geo.push_back(B[0]);
- geo.push_back(B[1]);
- geo.push_back(B[2]);
+ { // initialize field in the target region
+ assert(CbmKF::Instance()->vTargets.size() > 0);
+ auto& target = CbmKF::Instance()->vTargets[0];
+ algo->fCbmTargetX = target.x;
+ algo->fCbmTargetY = target.y;
+ algo->fCbmTargetZ = target.z;
+
+ for (int i = 0; i < 3; i++) {
+ Double_t point[3] = {0., 0., target.z + 2.5 * i};
+ Double_t B[3] = {0., 0., 0.};
+ if (CbmKF::Instance()->GetMagneticField()) { CbmKF::Instance()->GetMagneticField()->GetFieldValue(point, B); }
+ geo.push_back(point[2]);
+ geo.push_back(B[0]);
+ geo.push_back(B[1]);
+ geo.push_back(B[2]);
+ }
}
@@ -528,7 +541,7 @@ InitStatus CbmL1::Init()
z = t.z;
Xmax = Ymax = t.R;
- LOG(info) << "L1: Mvd station " << ist << " at z " << t.z << endl;
+ LOG(info) << "L1: Mvd station " << ist << " at z " << t.z;
}
@@ -562,7 +575,7 @@ InitStatus CbmL1::Init()
Xmax = station->GetXmax();
Ymax = station->GetYmax();
- LOG(info) << "L1: Sts station " << ist - NMvdStations << " at z " << station->GetZ() << endl;
+ LOG(info) << "L1: Sts station " << ist - NMvdStations << " at z " << station->GetZ();
}
if ((ist < (NMvdStations + NStsStations + NMuchStations)) && (ist >= (NMvdStations + NStsStations))) {
@@ -779,7 +792,7 @@ InitStatus CbmL1::Init()
const float RMax = hStaRadLen->GetXaxis()->GetXmax(); // should be same as min
algo->fRadThick[iSta].SetBins(NBins, RMax);
algo->fRadThick[iSta].table.resize(NBins);
-
+ double sumRL = 0., nRL = 0.;
for (int iB = 0; iB < NBins; iB++) {
algo->fRadThick[iSta].table[iB].resize(NBins);
for (int iB2 = 0; iB2 < NBins; iB2++) {
@@ -792,8 +805,11 @@ InitStatus CbmL1::Init()
// Correction for the incorrect harcoded value of RadThick of MVD stations
if (algo->fRadThick[iSta].table[iB][iB2] < 0.0015) algo->fRadThick[iSta].table[iB][iB2] = 0.0015;
// algo->fRadThick[iSta].table[iB][iB2] = algo->vStations[iSta].materialInfo.RadThick[0];
+ sumRL += algo->fRadThick[iSta].table[iB][iB2];
+ nRL++;
}
}
+ cout << " iSta " << iSta << " average RadLength in the material map " << sumRL / nRL << endl;
}
rlFile->Close();
rlFile->Delete();
@@ -831,7 +847,7 @@ InitStatus CbmL1::Init()
const float RMax = hStaRadLen->GetXaxis()->GetXmax(); // should be same as min
algo->fRadThick[iSta].SetBins(NBins, RMax);
algo->fRadThick[iSta].table.resize(NBins);
-
+ double sumRL = 0., nRL = 0.;
for (int iB = 0; iB < NBins; iB++) {
algo->fRadThick[iSta].table[iB].resize(NBins);
for (int iB2 = 0; iB2 < NBins; iB2++) {
@@ -839,8 +855,11 @@ InitStatus CbmL1::Init()
if (algo->fRadThick[iSta].table[iB][iB2] < algo->vStations[iSta].materialInfo.RadThick[0])
algo->fRadThick[iSta].table[iB][iB2] = algo->vStations[iSta].materialInfo.RadThick[0];
// cout << " iSta " << iSta << " iB " << iB << "iB2 " << iB2 << " RadThick " << algo->fRadThick[iSta].table[iB][iB2] << endl;
+ sumRL += algo->fRadThick[iSta].table[iB][iB2];
+ nRL++;
}
}
+ cout << " iSta " << iSta << " average RadLength in the material map " << sumRL / nRL << endl;
}
rlFile->Close();
rlFile->Delete();
diff --git a/reco/L1/CbmL1.h b/reco/L1/CbmL1.h
index baa9217323e43bcd0e2e245fa4f27e5ed415a74f..4428b32d51e4aac882c1198d9cd5f4bcb666c554 100644
--- a/reco/L1/CbmL1.h
+++ b/reco/L1/CbmL1.h
@@ -27,21 +27,17 @@
#include "CbmL1MCTrack.h"
#include "CbmL1Track.h"
#include "CbmL1Vtx.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"
-#include "L1Algo/L1Algo.h"
-#include "L1Algo/L1Vector.h"
-
-#include "CbmMCDataArray.h"
-#include "CbmMCEventList.h"
-#include "CbmTimeSlice.h"
-
#include "TClonesArray.h"
#include "TH1.h"
@@ -52,6 +48,8 @@
#include <set>
#include <utility>
+#include "L1Algo/L1Algo.h"
+#include "L1Algo/L1Vector.h"
#include "L1EventEfficiencies.h"
class L1AlgoInputData;
@@ -62,6 +60,7 @@ class L1FieldSlice;
class CbmL1Track;
class CbmL1MCTrack;
class KFTopoPerformance;
+class CbmMCDataObject;
class CbmEvent;
class CbmStsParSetSensor;
@@ -175,6 +174,7 @@ private:
// static bool compareZ(const int &a, const int &b );
// bool compareZ(const int &a, const int &b );
void Fill_vMCTracks();
+
/// Input Performance
void HitMatch(); // Procedure for match hits and MCPoints.
void FieldApproxCheck(); // Build histos with difference between Field map and approximated field
@@ -254,8 +254,6 @@ private:
Bool_t fUseTRD {false}; // if Trd data should be processed
Bool_t fUseTOF {false}; // if Tof data should be processed
- CbmL1Vtx PrimVtx {};
-
/// Input data
CbmTimeSlice* fTimeSlice {nullptr};
CbmMCEventList* fEventList {nullptr}; //! MC event list (all)
@@ -263,6 +261,7 @@ private:
CbmMCDataArray* fStsPoints {nullptr};
CbmMCDataArray* fMvdPoints {nullptr};
CbmMCDataArray* fMCTracks {nullptr};
+ CbmMCDataObject* fMcEventHeader {nullptr};
TClonesArray* listStsDigiMatch {nullptr};
TClonesArray* listStsClusters {nullptr};
diff --git a/reco/L1/CbmL1MCTrack.h b/reco/L1/CbmL1MCTrack.h
index 3e5867878672a619296ac7f29a9ef323e47d0d68..02458d19beef1f7df2e56c6c8fdd05fb28d431d8 100644
--- a/reco/L1/CbmL1MCTrack.h
+++ b/reco/L1/CbmL1MCTrack.h
@@ -84,6 +84,7 @@ public:
int iEvent = -1;
int mother_ID = -1;
int pdg = -1;
+ bool isSignal {0};
L1Vector<int> Points {"CbmL1MCTrack::Points"}; // indices of pints in L1::vMCPoints
L1Vector<int> StsHits {"CbmL1MCTrack::StsHits"}; // indices of hits in algo->vStsHits or L1::vStsHits
diff --git a/reco/L1/CbmL1Performance.cxx b/reco/L1/CbmL1Performance.cxx
index dee1d9b5a0d8d251c88d865d9a8544299396f800..df24b3afc0f15537398ad8fd0c574a89079e4285 100644
--- a/reco/L1/CbmL1Performance.cxx
+++ b/reco/L1/CbmL1Performance.cxx
@@ -369,7 +369,7 @@ void CbmL1::EfficienciesPerformance()
ntra.Inc(reco, killed, ratio_length, ratio_fakes, nclones, mc_length, mc_length_hits, "total");
- if ((mtra.IsPrimary()) && (mtra.z > 0)) { // D0
+ if (mtra.isSignal) { // D0
ntra.Inc(reco, killed, ratio_length, ratio_fakes, nclones, mc_length, mc_length_hits, "d0");
}
@@ -980,7 +980,7 @@ void CbmL1::HistoPerformance() // TODO: check if works correctly. Change vHitRe
}
- if ((mtra.IsPrimary()) && (mtra.z > 0)) { // D0
+ if (mtra.isSignal) { // D0
h_reco_d0_mom->Fill(momentum);
if (reco) p_eff_d0_vs_mom->Fill(momentum, 100.0);
else
diff --git a/reco/L1/CbmL1ReadEvent.cxx b/reco/L1/CbmL1ReadEvent.cxx
index 8ef2c32013d23700debfefda6bdcb890228d2a65..92aee4dbe902975739f33a7acec325c226e9c3aa 100644
--- a/reco/L1/CbmL1ReadEvent.cxx
+++ b/reco/L1/CbmL1ReadEvent.cxx
@@ -22,6 +22,7 @@
#include "CbmEvent.h"
#include "CbmKF.h"
#include "CbmL1.h"
+#include "CbmMCDataObject.h"
#include "CbmMatch.h"
#include "CbmMuchGeoScheme.h"
#include "CbmMuchPixelHit.h"
@@ -38,6 +39,8 @@
#include "CbmTrdHit.h"
#include "CbmTrdPoint.h"
+#include "FairMCEventHeader.h"
+
#include "L1Algo/L1Algo.h"
#include "L1AlgoInputData.h"
@@ -168,15 +171,22 @@ void CbmL1::ReadEvent(L1AlgoInputData* fData_, float& TsStart, float& TsLength,
if (fMvdPoints) {
Int_t nMvdPointsInEvent = fMvdPoints->Size(iFile, iEvent);
+ double maxDeviation = 0;
for (Int_t iMC = 0; iMC < nMvdPointsInEvent; iMC++) {
CbmL1MCPoint MC;
if (!ReadMCPoint(&MC, iMC, iFile, iEvent, 1)) {
- MC.iStation = -1;
- L1Station* sta = algo->vStations;
- for (Int_t iSt = 0; iSt < NStsStations; iSt++) {
- if (MC.z > sta[iSt].z[0] - 1) { MC.iStation = iSt; }
+ MC.iStation = -1;
+ L1Station* sta = algo->vStations;
+ double bestDist = 1.e20;
+ for (Int_t iSt = 0; iSt < NMvdStations; iSt++) {
+ double d = (MC.z - sta[iSt].z[0]);
+ if (fabs(d) < fabs(bestDist)) {
+ bestDist = d;
+ MC.iStation = iSt;
+ }
}
assert(MC.iStation >= 0);
+ if (fabs(maxDeviation) < fabs(bestDist)) { maxDeviation = bestDist; }
Double_t dtrck = dFEI(iFile, iEvent, MC.ID);
DFEI2I::iterator trk_it = dFEI2vMCTracks.find(dtrck);
assert(trk_it != dFEI2vMCTracks.end());
@@ -188,19 +198,29 @@ void CbmL1::ReadEvent(L1AlgoInputData* fData_, float& TsStart, float& TsLength,
nMvdPoints++;
}
}
+ // assert(fabs(maxDeviation)<1.);
+ if (fVerbose) { std::cout << "CbmL1ReadEvent: max deviation of Mvd points " << maxDeviation << std::endl; }
}
+
if (fStsPoints) {
- Int_t nMC = fStsPoints->Size(iFile, iEvent);
+ Int_t nMC = fStsPoints->Size(iFile, iEvent);
+ double maxDeviation = 0;
for (Int_t iMC = 0; iMC < nMC; iMC++) {
CbmL1MCPoint MC;
if (!ReadMCPoint(&MC, iMC, iFile, iEvent, 0)) {
- MC.iStation = -1;
- L1Station* sta = algo->vStations + NMvdStations;
+ MC.iStation = -1;
+ L1Station* sta = algo->vStations + NMvdStations;
+ double bestDist = 1.e20;
for (Int_t iSt = 0; iSt < NStsStations; iSt++) {
- if (MC.z > sta[iSt].z[0] - 2.5) { MC.iStation = NMvdStations + iSt; }
+ double d = (MC.z - sta[iSt].z[0]);
+ if (fabs(d) < fabs(bestDist)) {
+ bestDist = d;
+ MC.iStation = NMvdStations + iSt;
+ }
}
assert(MC.iStation >= 0);
+ if (fabs(maxDeviation) < fabs(bestDist)) { maxDeviation = bestDist; }
Double_t dtrck = dFEI(iFile, iEvent, MC.ID);
DFEI2I::iterator trk_it = dFEI2vMCTracks.find(dtrck);
assert(trk_it != dFEI2vMCTracks.end());
@@ -212,6 +232,8 @@ void CbmL1::ReadEvent(L1AlgoInputData* fData_, float& TsStart, float& TsLength,
nStsPoints++;
}
}
+ assert(fabs(maxDeviation) < 1.);
+ if (fVerbose) { std::cout << "CbmL1ReadEvent: max deviation of Sts points " << maxDeviation << std::endl; }
}
if (fMuchPoints) {
@@ -1064,81 +1086,64 @@ void CbmL1::ReadEvent(L1AlgoInputData* fData_, float& TsStart, float& TsLength,
void CbmL1::Fill_vMCTracks()
{
- PrimVtx.MC_ID = 999;
+ vMCTracks.clear();
{
- CbmL1Vtx Vtxcurr;
- int nvtracks = 0, nvtrackscurr = 0;
-
- vMCTracks.clear();
- {
- Int_t nMCTracks = 0;
- for (DFSET::iterator set_it = vFileEvent.begin(); set_it != vFileEvent.end(); ++set_it) {
- Int_t iFile = set_it->first;
- Int_t iEvent = set_it->second;
- nMCTracks += fMCTracks->Size(iFile, iEvent);
- }
- vMCTracks.reserve(nMCTracks);
- }
-
+ Int_t nMCTracks = 0;
for (DFSET::iterator set_it = vFileEvent.begin(); set_it != vFileEvent.end(); ++set_it) {
Int_t iFile = set_it->first;
Int_t iEvent = set_it->second;
+ nMCTracks += fMCTracks->Size(iFile, iEvent);
+ }
+ vMCTracks.reserve(nMCTracks);
+ }
- Int_t nMCTrack = fMCTracks->Size(iFile, iEvent);
+ int fileEvent = 0;
+ for (DFSET::iterator set_it = vFileEvent.begin(); set_it != vFileEvent.end(); ++set_it, ++fileEvent) {
+ Int_t iFile = set_it->first;
+ Int_t iEvent = set_it->second;
- for (Int_t iMCTrack = 0; iMCTrack < nMCTrack; iMCTrack++) {
- CbmMCTrack* MCTrack = L1_DYNAMIC_CAST<CbmMCTrack*>(fMCTracks->Get(iFile, iEvent, iMCTrack));
- if (!MCTrack) continue;
+ auto header = dynamic_cast<FairMCEventHeader*>(fMcEventHeader->Get(iFile, iEvent));
+ assert(header);
+ if (fVerbose) {
+ cout << "mc event vertex at " << header->GetX() << " " << header->GetY() << " " << header->GetZ() << endl;
+ }
- int mother_ID = MCTrack->GetMotherId();
+ Int_t nMCTrack = fMCTracks->Size(iFile, iEvent);
+ for (Int_t iMCTrack = 0; iMCTrack < nMCTrack; iMCTrack++) {
+ CbmMCTrack* MCTrack = L1_DYNAMIC_CAST<CbmMCTrack*>(fMCTracks->Get(iFile, iEvent, iMCTrack));
+ if (!MCTrack) continue;
- if (mother_ID < 0 && mother_ID != -2) mother_ID = -iEvent - 1;
- TVector3 vr;
- TLorentzVector vp;
- MCTrack->GetStartVertex(vr);
- MCTrack->Get4Momentum(vp);
+ int mother_ID = MCTrack->GetMotherId();
- Int_t pdg = MCTrack->GetPdgCode();
- Double_t q = 1, mass = 0.;
- if (pdg < 9999999 && ((TParticlePDG*) TDatabasePDG::Instance()->GetParticle(pdg))) {
- q = TDatabasePDG::Instance()->GetParticle(pdg)->Charge() / 3.0;
- mass = TDatabasePDG::Instance()->GetParticle(pdg)->Mass();
- }
- else
- q = 0;
- Int_t IND_Track = vMCTracks.size(); //or iMCTrack?
- CbmL1MCTrack T(mass, q, vr, vp, IND_Track, mother_ID, pdg);
- // CbmL1MCTrack T(mass, q, vr, vp, iMCTrack, mother_ID, pdg);
- T.time = MCTrack->GetStartT();
- T.iFile = iFile;
- T.iEvent = iEvent;
-
- vMCTracks.push_back(T);
- // Double_t dtrck =dFEI(iFile,iEvent,iMCTrack);
- dFEI2vMCTracks.insert(DFEI2I::value_type(dFEI(iFile, iEvent, iMCTrack), vMCTracks.size() - 1));
-
- if (T.mother_ID < 0) { // vertex track
- if (PrimVtx.MC_ID == 999 || fabs(T.z - Vtxcurr.MC_z) > 1.e-7) { // new vertex
- if (nvtrackscurr > nvtracks) {
- PrimVtx = Vtxcurr;
- nvtracks = nvtrackscurr;
- }
- Vtxcurr.MC_x = T.x;
- Vtxcurr.MC_y = T.y;
- Vtxcurr.MC_z = T.z;
- Vtxcurr.MC_ID = T.mother_ID;
- nvtrackscurr = 1;
- }
- else
- nvtrackscurr++;
- }
- } //iTracks
- } //Links
- if (nvtrackscurr > nvtracks) PrimVtx = Vtxcurr;
- } //PrimVtx
+ if (mother_ID < 0 && mother_ID != -2) mother_ID = -iEvent - 1;
+ TVector3 vr;
+ TLorentzVector vp;
+ MCTrack->GetStartVertex(vr);
+ MCTrack->Get4Momentum(vp);
+
+ Int_t pdg = MCTrack->GetPdgCode();
+ Double_t q = 0, mass = 0.;
+ if (pdg < 9999999 && ((TParticlePDG*) TDatabasePDG::Instance()->GetParticle(pdg))) {
+ q = TDatabasePDG::Instance()->GetParticle(pdg)->Charge() / 3.0;
+ mass = TDatabasePDG::Instance()->GetParticle(pdg)->Mass();
+ }
+
+ Int_t IND_Track = vMCTracks.size(); //or iMCTrack?
+ CbmL1MCTrack T(mass, q, vr, vp, IND_Track, mother_ID, pdg);
+ // CbmL1MCTrack T(mass, q, vr, vp, iMCTrack, mother_ID, pdg);
+ T.time = MCTrack->GetStartT();
+ T.iFile = iFile;
+ T.iEvent = iEvent;
+ // signal: primary tracks, displaced from the primary vertex
+ T.isSignal = T.IsPrimary() && (T.z > header->GetZ() + 1.e-10);
+
+ vMCTracks.push_back(T);
+ // Double_t dtrck =dFEI(iFile,iEvent,iMCTrack);
+ dFEI2vMCTracks.insert(DFEI2I::value_type(dFEI(iFile, iEvent, iMCTrack), vMCTracks.size() - 1));
+ } //iTracks
+ } //Links
- if (fVerbose && PrimVtx.MC_ID == 999) cout << "No primary vertex !!!" << endl;
} //Fill_vMCTracks
bool CbmL1::ReadMCPoint(CbmL1MCPoint* MC, int iPoint, int file, int event, int MVD)
diff --git a/reco/L1/L1Algo/L1Algo.cxx b/reco/L1/L1Algo/L1Algo.cxx
index b806bfd33995fdbec14bdce88129804c48e94812..3c2c4f3c151298eaf4a63dfdf8d16e82f86968c2 100644
--- a/reco/L1/L1Algo/L1Algo.cxx
+++ b/reco/L1/L1Algo/L1Algo.cxx
@@ -296,8 +296,8 @@ void L1Algo::GetHitCoor(const L1Hit& _h, fscal& _x, fscal& _y, char iS)
L1Station& sta = vStations[int(iS)];
fscal u = _h.u;
fscal v = _h.v;
- _x = (sta.xInfo.sin_phi[0] * u + sta.xInfo.cos_phi[0] * v) / _h.z;
- _y = (sta.yInfo.cos_phi[0] * u + sta.yInfo.sin_phi[0] * v) / _h.z;
+ _x = (sta.xInfo.sin_phi[0] * u + sta.xInfo.cos_phi[0] * v) / (_h.z - fCbmTargetZ[0]);
+ _y = (sta.yInfo.cos_phi[0] * u + sta.yInfo.sin_phi[0] * v) / (_h.z - fCbmTargetZ[0]);
}
void L1Algo::GetHitCoor(const L1Hit& _h, fscal& _x, fscal& _y, fscal& _z, const L1Station& sta)
diff --git a/reco/L1/L1Algo/L1Algo.h b/reco/L1/L1Algo/L1Algo.h
index b445714872c96af069657b087ce38b193ccb47e0..7b4f244a78554358684b28aeed5e9c1a560545b9 100644
--- a/reco/L1/L1Algo/L1Algo.h
+++ b/reco/L1/L1Algo/L1Algo.h
@@ -39,6 +39,7 @@ class L1AlgoDraw;
#include <iomanip>
#include <iostream>
+#include <limits>
#include <map>
#include "L1Branch.h"
@@ -627,6 +628,7 @@ private:
map<int, int> threadNumberToCpuMap {};
+ static constexpr float kNaN {std::numeric_limits<float>::signaling_NaN()};
float TRACK_CHI2_CUT {10.f};
float TRIPLET_CHI2_CUT {5.f}; // cut for selecting triplets before collecting tracks.per one DoF
@@ -635,18 +637,22 @@ private:
float TIME_CUT2 {0.f};
fvec MaxDZ {
- 0.f}; // correction in order to take into account overlaping and iff z. if sort by y then it is max diff between same station's modules (~0.4cm)
+ kNaN}; // correction in order to take into account overlaping and iff z. if sort by y then it is max diff between same station's modules (~0.4cm)
/// parameters which are different for different iterations. Set in the begin of CAL1TrackFinder
- float Pick_gather {0.f}; // same for attaching additional hits to track
- float PickNeighbour {0.f}; // (PickNeighbour < dp/dp_error) => triplets are neighbours
- fvec MaxInvMom {0.f}; // max considered q/p for tracks
- fvec MaxSlopePV {0.f}; // max slope (tx\ty) in prim vertex
- float MaxSlope {0.f}; // max slope (tx\ty) in 3d hit position of a triplet
- fvec targX {0.f}; // target coor TODO: set defaults to a crasy value to be sure that the target is initialised later
- fvec targY {0.f};
- fvec targZ {0.f};
+ float Pick_gather {kNaN}; // same for attaching additional hits to track
+ float PickNeighbour {kNaN}; // (PickNeighbour < dp/dp_error) => triplets are neighbours
+ fvec MaxInvMom {kNaN}; // max considered q/p for tracks
+ fvec MaxSlopePV {kNaN}; // max slope (tx\ty) in prim vertex
+ float MaxSlope {kNaN}; // max slope (tx\ty) in 3d hit position of a triplet
+ fvec fCbmTargetX {kNaN}; // target position
+ fvec fCbmTargetY {kNaN};
+ fvec fCbmTargetZ {kNaN};
+ fvec fTargX {kNaN}; // target position for the current iteration
+ fvec fTargY {kNaN};
+ fvec fTargZ {kNaN};
+
L1FieldValue targB _fvecalignment {}; // field in the target point
L1XYMeasurementInfo TargetXYInfo _fvecalignment {}; // target constraint [cm]
diff --git a/reco/L1/L1Algo/L1CATrackFinder.cxx b/reco/L1/L1Algo/L1CATrackFinder.cxx
index 4c1275bb8b3a3074401f6da5ed4171c64802393f..e3d5e11b28fc42f187e15ae0f05a29b109aead5e 100644
--- a/reco/L1/L1Algo/L1CATrackFinder.cxx
+++ b/reco/L1/L1Algo/L1CATrackFinder.cxx
@@ -132,7 +132,7 @@ inline void L1Algo::f11( /// input 1st stage of singlet search
fvec zstal = stal.z;
fvec zstam = stam.z;
- int istal_global = 5;
+ int istal_global = 5; //TODO: SG: what are these stations? Should the numbers depend on the presence of mvd?
int istam_global = 6;
L1Station& stal_global = vStations[istal_global];
L1Station& stam_global = vStations[istam_global];
@@ -168,8 +168,7 @@ inline void L1Algo::f11( /// input 1st stage of singlet search
fvec zl = zPos_l[i1_V];
fvec& time = HitTime_l[i1_V];
fvec& timeEr = HitTimeEr[i1_V];
- const fvec dzli = 1. / (zl - targZ);
-
+ const fvec dzli = 1.f / (zl - fTargZ);
fvec dx1, dy1, dxy1 = 0;
dUdV_to_dX(d_u[i1_V], d_v[i1_V], dx1, stal);
@@ -178,8 +177,8 @@ inline void L1Algo::f11( /// input 1st stage of singlet search
StripsToCoor(u, v, xl, yl, stal);
- const fvec tx = (xl - targX) * dzli;
- const fvec ty = (yl - targY) * dzli;
+ const fvec tx = (xl - fTargX) * dzli;
+ const fvec ty = (yl - fTargY) * dzli;
// estimate field for singlet creation
int istac = istal / 2; // "center" station
@@ -193,31 +192,31 @@ inline void L1Algo::f11( /// input 1st stage of singlet search
L1Station& stac_global = vStations[istac_global];
fvec zstac_global = stac.z;
- stac.fieldSlice.GetFieldValue(targX + tx * (zstac - targZ), targY + ty * (zstac - targZ), centB);
- stal.fieldSlice.GetFieldValue(targX + tx * (zstal - targZ), targY + ty * (zstal - targZ), l_B);
+ stac.fieldSlice.GetFieldValue(fTargX + tx * (zstac - fTargZ), fTargY + ty * (zstac - fTargZ), centB);
+ stal.fieldSlice.GetFieldValue(fTargX + tx * (zstal - fTargZ), fTargY + ty * (zstal - fTargZ), l_B);
- stam_global.fieldSlice.GetFieldValue(targX + tx * (zstam_global - targZ), targY + ty * (zstam_global - targZ),
+ stam_global.fieldSlice.GetFieldValue(fTargX + tx * (zstam_global - fTargZ), fTargY + ty * (zstam_global - fTargZ),
m_B_global);
- stal_global.fieldSlice.GetFieldValue(targX + tx * (zstal_global - targZ), targY + ty * (zstal_global - targZ),
+ stal_global.fieldSlice.GetFieldValue(fTargX + tx * (zstal_global - fTargZ), fTargY + ty * (zstal_global - fTargZ),
l_B_global);
- stac_global.fieldSlice.GetFieldValue(targX + tx * (zstac_global - targZ), targY + ty * (zstac_global - targZ),
+ stac_global.fieldSlice.GetFieldValue(fTargX + tx * (zstac_global - fTargZ), fTargY + ty * (zstac_global - fTargZ),
centB_global);
- if (istac != istal) fld0.Set(l_B, stal.z, centB, stac.z, targB, targZ);
+ if (istac != istal) fld0.Set(l_B, stal.z, centB, stac.z, targB, fTargZ);
else
- fld0.Set(l_B, zstal, targB, targZ);
+ fld0.Set(l_B, zstal, targB, fTargZ);
// estimate field for the next extrapolations
- stam.fieldSlice.GetFieldValue(targX + tx * (zstam - targZ), targY + ty * (zstam - targZ), m_B);
+ stam.fieldSlice.GetFieldValue(fTargX + tx * (zstam - fTargZ), fTargY + ty * (zstam - fTargZ), m_B);
#define USE_3HITS
#ifndef USE_3HITS
if (istac != istal) fld1.Set(m_B, zstam, l_B, zstal, centB, zstac);
else
- fld1.Set(m_B, zstam, l_B, zstal, targB, targZ);
+ fld1.Set(m_B, zstam, l_B, zstal, targB, fTargZ);
#else // if USE_3HITS // the best now
L1FieldValue r_B _fvecalignment;
L1Station& star = vStations[istam + 1];
fvec zstar = star.z;
- star.fieldSlice.GetFieldValue(targX + tx * (zstar - targZ), targY + ty * (zstar - targZ), r_B);
+ star.fieldSlice.GetFieldValue(fTargX + tx * (zstar - fTargZ), fTargY + ty * (zstar - fTargZ), r_B);
fld1.Set(r_B, zstar, m_B, zstam, l_B, zstal);
if ((istam + 1) >= fNfieldStations)
fld1.Set(m_B_global, zstam_global, l_B_global, zstal_global, centB_global, zstac_global);
@@ -265,7 +264,7 @@ inline void L1Algo::f11( /// input 1st stage of singlet search
if (istal < fNfieldStations) {
fvec eX, eY, J04, J14;
fvec dz;
- dz = targZ - zl;
+ dz = fTargZ - zl;
L1ExtrapolateJXY0(T.tx, T.ty, dz, fld0, eX, eY, J04, J14);
fvec J[6];
J[0] = dz;
@@ -274,11 +273,11 @@ inline void L1Algo::f11( /// input 1st stage of singlet search
J[3] = 0;
J[4] = dz;
J[5] = J14;
- L1FilterVtx(T, targX, targY, TargetXYInfo, eX, eY, J);
+ L1FilterVtx(T, fTargX, fTargY, TargetXYInfo, eX, eY, J);
}
else {
- L1ExtrapolateLine(T, targZ);
- L1FilterXY(T, targX, targY, TargetXYInfo);
+ L1ExtrapolateLine(T, fTargZ);
+ L1FilterXY(T, fTargX, fTargY, TargetXYInfo);
}
}
@@ -289,7 +288,7 @@ inline void L1Algo::f11( /// input 1st stage of singlet search
{
fvec eX, eY, J04, J14;
fvec dz;
- dz = targZ - zl;
+ dz = fTargZ - zl;
L1ExtrapolateJXY0(T.tx, T.ty, dz, fld0, eX, eY, J04, J14);
fvec J[6];
J[0] = dz;
@@ -298,16 +297,16 @@ inline void L1Algo::f11( /// input 1st stage of singlet search
J[3] = 0;
J[4] = dz;
J[5] = J14;
- L1FilterVtx(T, targX, targY, TargetXYInfo, eX, eY, J);
+ L1FilterVtx(T, fTargX, fTargY, TargetXYInfo, eX, eY, J);
}
#else // TODO: doesn't work. Why?
- T.x = targX;
- T.y = targY;
+ T.x = fTargX;
+ T.y = fTargY;
- T.z = targZ;
+ T.z = fTargZ;
T.C00 = TargetXYInfo.C00;
T.C10 = TargetXYInfo.C10;
T.C11 = TargetXYInfo.C11;
@@ -411,7 +410,7 @@ inline void L1Algo::f20(
// Pick_m22 is not used, search for mean squared, 2nd version
// -- collect possible doublets --
- const fscal iz = 1 / T1.z[i1_4];
+ const fscal iz = 1.f / (T1.z[i1_4] - fCbmTargetZ[0]);
const float& timeError = T1.C55[i1_4];
const float& time = T1.t[i1_4];
@@ -719,7 +718,7 @@ inline void L1Algo::f30( // input
#ifdef DO_NOT_SELECT_TRIPLETS
if (isec == TRACKS_FROM_TRIPLETS_ITERATION) Pick_r22 = Pick_r2 + 1;
#endif
- const fscal iz = 1 / T2.z[i2_4];
+ const fscal iz = 1.f / (T2.z[i2_4] - fCbmTargetZ[0]);
L1HitAreaTime area(vGridTime[&star - vStations], T2.x[i2_4] * iz, T2.y[i2_4] * iz,
(sqrt(Pick_r22 * (T2.C00 + stam.XYInfo.C00)) + MaxDZ * fabs(T2.tx))[i2_4] * iz,
(sqrt(Pick_r22 * (T2.C11 + stam.XYInfo.C11)) + MaxDZ * fabs(T2.ty))[i2_4] * iz, time,
@@ -1860,9 +1859,9 @@ void L1Algo::CATrackFinder()
MaxSlope = 2.748; // corresponds to 70 grad
// define the target
- targX = 0;
- targY = 0;
- targZ = 0; // suppose, what target will be at (0,0,0)
+ fTargX = fCbmTargetX;
+ fTargY = fCbmTargetY;
+ fTargZ = fCbmTargetZ;
float SigmaTargetX = 0, SigmaTargetY = 0; // target constraint [cm]
if ((isec == kFastPrimIter) || (isec == kFastPrimIter2) || (isec == kFastPrimJumpIter) || (isec == kAllPrimIter)
@@ -1877,7 +1876,7 @@ void L1Algo::CATrackFinder()
|| (isec == kAllSecJumpIter)) { //use outer radius of the 1st station as a constraint
L1Station& st = vStations[0];
SigmaTargetX = SigmaTargetY = 10; //st.Rmax[0];
- targZ = 0.; //-1.;
+ fTargZ = fCbmTargetZ; // fCbmTargetZ-1.;
st.fieldSlice.GetFieldValue(0, 0, targB);
}
@@ -2358,7 +2357,7 @@ void L1Algo::CATrackFinder()
float xcoor, ycoor = 0;
L1Station stah = vStations[0];
StripsToCoor(tempPoint.U(), tempPoint.V(), xcoor, ycoor, stah);
- float zcoor = tempPoint.Z();
+ float zcoor = tempPoint.Z() - fCbmTargetZ[0];
float timeFlight = sqrt(xcoor * xcoor + ycoor * ycoor + zcoor * zcoor) / 30.f; // c = 30[cm/ns]
sumTime += (hit.t - timeFlight);
diff --git a/reco/L1/L1Algo/L1HitPoint.h b/reco/L1/L1Algo/L1HitPoint.h
index c03eb6d67f1aeba4ae25b5bd96c8b1a57c722258..ceb27e5c918a7cf2826e2ac525471718a4ba82ca 100644
--- a/reco/L1/L1Algo/L1HitPoint.h
+++ b/reco/L1/L1Algo/L1HitPoint.h
@@ -99,6 +99,7 @@ public:
static const float R = 60;
static const float shortPackingConstant = 2 * R / 65535.f;
+//TODO: change the Z conversion. Hit Z maybe negative in the new setup.
static const float MZ = 110;
static const float shortPackingConstantZ = MZ / 65535.f;
diff --git a/reco/L1/L1Algo/L1Track.h b/reco/L1/L1Algo/L1Track.h
index 181b6b3aeba866c80f72a6da57ccd4797c64fb05..1cd1ea6216d89077a882ea7b05f67805c3cf8d3d 100644
--- a/reco/L1/L1Algo/L1Track.h
+++ b/reco/L1/L1Algo/L1Track.h
@@ -21,12 +21,17 @@
#ifndef L1Track_H
#define L1Track_H
+#include <limits>
+
class L1Track {
public:
+ static constexpr float kNaN {std::numeric_limits<float>::signaling_NaN()};
+
unsigned char NHits;
unsigned char n;
- float Momentum, fTrackTime;
- fscal TFirst[7], CFirst[21], TLast[7], CLast[21], Tpv[7], Cpv[21], chi2;
+ float Momentum {kNaN}, fTrackTime {kNaN};
+ fscal TFirst[7] {kNaN}, CFirst[21] {kNaN}, TLast[7] {kNaN}, CLast[21] {kNaN}, Tpv[7] {kNaN}, Cpv[21] {kNaN},
+ chi2 {kNaN};
short int NDF;
int FirstHitIndex, LastHitIndex;
diff --git a/reco/detectors/sts/qa/CbmStsFindTracksQa.cxx b/reco/detectors/sts/qa/CbmStsFindTracksQa.cxx
index b60ea95b0293ea8a65be6ef04c124ef9ed252513..ff51d163886e13059f755fc8061463b32767c0e2 100644
--- a/reco/detectors/sts/qa/CbmStsFindTracksQa.cxx
+++ b/reco/detectors/sts/qa/CbmStsFindTracksQa.cxx
@@ -491,7 +491,10 @@ InitStatus CbmStsFindTracksQa::GetGeometry()
// Get target geometry
GetTargetPosition();
fMvdNstations = 0;
- {
+ // CbmMvdDetector may issue an error when there is no MVD in the setup
+ // For now we determine the presence of MVD from the presence of the mvd hit branch
+ // TODO: use CbmSetup to determine if MVD is present
+ if (fManager->GetObject("MvdHit")) {
CbmMvdDetector* mvdDetector = CbmMvdDetector::Instance();
if (mvdDetector) {
CbmMvdStationPar* mvdStationPar = mvdDetector->GetParameterFile();