diff --git a/reco/alignment/CbmBbaAlignmentTask.cxx b/reco/alignment/CbmBbaAlignmentTask.cxx
index 536a380bbcb92c56337d8e4c70850ac79ecb8a02..3e854451a176d63b5e4f531154bb347126935d74 100644
--- a/reco/alignment/CbmBbaAlignmentTask.cxx
+++ b/reco/alignment/CbmBbaAlignmentTask.cxx
@@ -16,6 +16,7 @@
// ROOT headers
#include "FairRootManager.h"
+#include "FairRunAna.h"
#include "TClonesArray.h"
#include "TFile.h"
#include "TGeoPhysicalNode.h"
@@ -296,6 +297,8 @@ void CbmBbaAlignmentTask::Exec(Option_t* /*opt*/)
ca::Framework* l1 = CbmL1::Instance()->fpAlgo;
const ca::Parameters<ca::fvec>& l1Par = l1->GetParameters();
+ fFitter.SetDoSmooth(true);
+
// select tracks for alignment and store them
if (fTrackingMode == kSts && fInputStsTracks) {
@@ -334,7 +337,6 @@ void CbmBbaAlignmentTask::Exec(Option_t* /*opt*/)
fFitter.SetDefaultMomentumForMs(0.5);
fFitter.FixMomentumForMs(true);
- fFitter.SetDoSmooth(true);
for (int iTr = 0; iTr < fInputGlobalTracks->GetEntriesFast(); iTr++) {
if (static_cast<int>(fTracks.size()) >= fMaxNtracks) {
@@ -372,63 +374,60 @@ void CbmBbaAlignmentTask::Exec(Option_t* /*opt*/)
LOG(warning) << "failed to fit the track! ";
continue;
}
+ t.fAlignedTrack = t.fUnalignedTrack;
+ fTracks.push_back(t);
+ }
+ } // end of mcbm tracking mode
- // fix the material radiation thickness to avoid the chi2 rattling at the material map bin boundaries
- // (to be improved)
- {
- bool isGood = true;
- for (unsigned int in = 0; in < t.fUnalignedTrack.fNodes.size(); in++) {
- CbmKfTrackFitter::FitNode& node = t.fUnalignedTrack.fNodes[in];
- if (!node.fIsFitted) {
- isGood = false;
- break;
- }
- node.fIsRadThickFixed = true;
- if (node.fRadThick > 0.01) {
- LOG(warning) << "CbmBbaAlignmentTask: radiation thickness is too large: " << node.fRadThick;
- node.fRadThick = 0.01;
- }
- if (node.fRadThick < 0.001) {
- LOG(warning) << "CbmBbaAlignmentTask: radiation thickness is too small: " << node.fRadThick;
- node.fRadThick = 0.001;
- }
- }
- if (!isGood) continue;
- }
- // ensure that all the hits are not too much deviated from the trajectory
- // (to be improved)
- {
- bool isGood = true;
- const double cutX = 8.;
- const double cutY = 8.;
- for (unsigned int in = 0; in < t.fUnalignedTrack.fNodes.size(); in++) {
- CbmKfTrackFitter::Track tr = t.fUnalignedTrack;
- CbmKfTrackFitter::FitNode& node = tr.fNodes[in];
- if (!node.fIsXySet) {
- continue;
- }
- node.fIsXySet = false; // exclude the node from the fit
- if (!fFitter.FitTrack(tr)) {
- isGood = false;
- break;
- }
- double x_residual = node.fMxy.X() - node.fTrack.X(); // this should be the difference vector
- double y_residual = node.fMxy.Y() - node.fTrack.Y();
- double x_pull = x_residual / sqrt(node.fMxy.Dx2() + node.fTrack.GetCovariance(0, 0));
- double y_pull = y_residual / sqrt(node.fMxy.Dy2() + node.fTrack.GetCovariance(1, 1));
- if (!node.fIsFitted || (node.fMxy.NdfX() > 0. && fabs(x_pull) > cutX)
- || (node.fMxy.NdfY() > 0. && fabs(y_pull) > cutY)) {
- isGood = false;
- break;
- }
- }
- if (!isGood) continue;
+ // fix the material radiation thickness to avoid the chi2 rattling at the material map bin boundaries
+ // (to be improved)
+ for (TrackContainer& t : fTracks) {
+ if (!t.fIsActive) continue;
+ for (unsigned int in = 0; in < t.fUnalignedTrack.fNodes.size(); in++) {
+ CbmKfTrackFitter::FitNode& node = t.fUnalignedTrack.fNodes[in];
+ if (!node.fIsFitted) {
+ t.fIsActive = false;
+ break;
}
+ node.fIsRadThickFixed = true;
+ if (node.fRadThick > 0.01) {
+ LOG(warning) << "CbmBbaAlignmentTask: radiation thickness is too large: " << node.fRadThick;
+ node.fRadThick = 0.01;
+ }
+ if (node.fRadThick < 0.001) {
+ LOG(warning) << "CbmBbaAlignmentTask: radiation thickness is too small: " << node.fRadThick;
+ node.fRadThick = 0.001;
+ }
+ }
+ }
- t.fAlignedTrack = t.fUnalignedTrack;
-
- fTracks.push_back(t);
+ // ensure that all the hits are not too much deviated from the trajectory
+ // (to be improved)
+ for (TrackContainer& t : fTracks) {
+ if (!t.fIsActive) continue;
+ const double cutX = 8.;
+ const double cutY = 8.;
+ for (unsigned int in = 0; in < t.fUnalignedTrack.fNodes.size(); in++) {
+ CbmKfTrackFitter::Track tr = t.fUnalignedTrack;
+ CbmKfTrackFitter::FitNode& node = tr.fNodes[in];
+ if (!node.fIsXySet) {
+ continue;
+ }
+ node.fIsXySet = false; // exclude the node from the fit
+ if (!fFitter.FitTrack(tr)) {
+ t.fIsActive = false;
+ break;
+ }
+ double x_residual = node.fMxy.X() - node.fTrack.X(); // this should be the difference vector
+ double y_residual = node.fMxy.Y() - node.fTrack.Y();
+ double x_pull = x_residual / sqrt(node.fMxy.Dx2() + node.fTrack.GetCovariance(0, 0));
+ double y_pull = y_residual / sqrt(node.fMxy.Dy2() + node.fTrack.GetCovariance(1, 1));
+ if (!node.fIsFitted || (node.fMxy.NdfX() > 0. && fabs(x_pull) > cutX)
+ || (node.fMxy.NdfY() > 0. && fabs(y_pull) > cutY)) {
+ t.fIsActive = false;
+ break;
+ }
}
}
@@ -505,9 +504,15 @@ void CbmBbaAlignmentTask::ApplyAlignment(const std::vector<double>& par)
double dy = parConstrained[3 * s.fAlignmentBody + 1];
double dz = parConstrained[3 * s.fAlignmentBody + 2];
+ // shift the hit
nodeAligned.fMxy.SetX(node.fMxy.X() + dx);
nodeAligned.fMxy.SetY(node.fMxy.Y() + dy);
nodeAligned.fZ = node.fZ + dz;
+
+ // shift the fitted track to the Z position of the hit
+ nodeAligned.fTrack.SetX(node.fTrack.X() + node.fTrack.Tx() * dz);
+ nodeAligned.fTrack.SetY(node.fTrack.Y() + node.fTrack.Ty() * dz);
+ nodeAligned.fTrack.SetZ(node.fTrack.Z() + dz);
}
}
@@ -538,19 +543,29 @@ double CbmBbaAlignmentTask::CostFunction(const std::vector<double>& par)
std::vector<double> chi2Thread(fNthreads, 0.);
std::vector<long> ndfThread(fNthreads, 0);
- fFitter.SetDoSmooth(false);
-
auto fitThread = [&](int iThread) {
CbmKfTrackFitter fit(fFitter);
+ fit.SetDoSmooth(false);
+ fit.SetNoMultipleScattering();
+
for (unsigned int itr = iThread; itr < fTracks.size(); itr += fNthreads) {
auto& t = fTracks[itr];
if (!t.fIsActive) continue;
+
+ // fit.SetDebugInfo(" track " + std::to_string(itr));
+
+ for (unsigned int in = 0; in < t.fAlignedTrack.fNodes.size(); in++) {
+ CbmKfTrackFitter::FitNode& node = t.fAlignedTrack.fNodes[in];
+ if (!node.fIsFitted) {
+ LOG(fatal) << "BBA: Cost function: aligned node is not fitted! ";
+ }
+ }
bool ok = fit.FitTrack(t.fAlignedTrack);
double chi2 = t.fAlignedTrack.fNodes.back().fTrack.GetChiSq();
double ndf = t.fAlignedTrack.fNodes.back().fTrack.GetNdf();
if (!ok || !std::isfinite(chi2) || chi2 <= 0. || ndf <= 0.) {
+ // TODO: deactivate the track and continue
LOG(fatal) << "BBA: fit failed! ";
- assert(0);
chi2 = 0.;
ndf = 0.;
}
@@ -571,7 +586,6 @@ double CbmBbaAlignmentTask::CostFunction(const std::vector<double>& par)
th.join();
}
-
fChi2Total = 0.;
fNdfTotal = 0;
for (int i = 0; i < fNthreads; i++) {
@@ -663,7 +677,7 @@ void CbmBbaAlignmentTask::Finish()
t.fAlignedTrack = t.fUnalignedTrack;
}
- if (0) { // one alignment body per tracking station
+ if (1) { // one alignment body per tracking station
fNalignmentBodies = fNtrackingStations;
fAlignmentBodies.resize(fNalignmentBodies);
@@ -832,19 +846,26 @@ void CbmBbaAlignmentTask::Finish()
fCostInitial = CostFunction(parMisaligned);
fFixedNdf = -1; //fNdfTotal;
+ // plot the residuals before alignment
- for (auto& t : fTracks) {
+ for (const auto& t : fTracks) {
if (!t.fIsActive) continue;
// calculate the residuals for all tracks at each fitNode before alignment
for (unsigned int in = 0; in < t.fAlignedTrack.fNodes.size(); in++) {
CbmKfTrackFitter::Track tr = t.fAlignedTrack;
CbmKfTrackFitter::FitNode& node = tr.fNodes[in];
+
if (!node.fIsXySet) {
continue;
}
node.fIsXySet = false;
+
fFitter.FitTrack(tr);
+ if (!node.fIsFitted) {
+ continue;
+ }
+
Sensor& s = fSensors[node.fReference1];
double x_residual = node.fMxy.X() - node.fTrack.X(); // this should be the difference vector
@@ -877,20 +898,24 @@ void CbmBbaAlignmentTask::Finish()
LOG(info) << " Cost function for the true parameters: " << fCostIdeal;
LOG(info) << " Cost function for the initial parameters: " << fCostInitial;
+ // run the alignment
alignment.align();
double costResult = CostFunction(alignment.getResult());
-
CostFunction(alignment.getResult());
CostFunction(alignment.getResult());
// Create alignment matrices:
std::vector<double> result = alignment.getResult();
std::map<std::string, TGeoHMatrix> alignmentMatrices;
- for (int i = 0; i < fNalignmentBodies; i++) {
+
+ for (auto& s : fSensors) {
+ int i = s.fAlignmentBody;
+ assert(i < fNalignmentBodies);
+ if (i < 0) continue;
+
// For STS sensors
- auto& s = fSensors[i];
if (s.fSystemId == ECbmModuleId::kSts) {
// get sensor for sensorID
const CbmStsElement* el = CbmStsSetup::Instance()->GetElement(s.fSensorId, kStsSensor);
@@ -925,7 +950,7 @@ void CbmBbaAlignmentTask::Finish()
// alignmentMatrices.insert(AlignNode(Form("fTrackingStation %i", s.fTrackingStation), result[i * 3 + 0],
// result[i * 3 + 1], result[i * 3 + 2], 0., 0., 0.));
// }
- }
+ } // sensors
// save matrices to disk
TFile* misalignmentMatrixRootfile =