diff --git a/macro/alignment/run_BbaAlignment.C b/macro/alignment/run_BbaAlignment.C
index 4f1db8c479af950a9444cc25f849558333380b76..c8dad80b6f36027f4c9390e91c8f287408ba893f 100644
--- a/macro/alignment/run_BbaAlignment.C
+++ b/macro/alignment/run_BbaAlignment.C
@@ -43,7 +43,7 @@
#endif
void run_BbaAlignment(TString input = "", TString output = "", TString setup = "sis100_electron",
- TString paramFile = "", Int_t nEvents = -1)
+ TString paramFile = "", Int_t nEvents = -1, double misalignmentRange = 0., int seed = 1)
{
gROOT->SetBatch(kTRUE);
@@ -151,6 +151,8 @@ void run_BbaAlignment(TString input = "", TString output = "", TString setup = "
// ----- BBA alignment --------------------------------------------
CbmBbaAlignmentTask* alignment = new CbmBbaAlignmentTask();
alignment->SetStsTrackingMode();
+ alignment->SetSimulatedMisalignmentRange(misalignmentRange);
+ alignment->SetRandomSeed(seed);
run->AddTask(alignment);
diff --git a/reco/alignment/CbmBbaAlignmentTask.cxx b/reco/alignment/CbmBbaAlignmentTask.cxx
index 1f4aa751ba48aef297138aad399f05fd77b067b6..255b7686d01e68208da66026bf53e8622a011767 100644
--- a/reco/alignment/CbmBbaAlignmentTask.cxx
+++ b/reco/alignment/CbmBbaAlignmentTask.cxx
@@ -92,6 +92,11 @@ void CbmBbaAlignmentTask::TrackContainer::MakeConsistent()
}
}
+// void addVirtualMisalignment()
+// {
+
+
+// }
CbmBbaAlignmentTask::CbmBbaAlignmentTask(const char* name, Int_t iVerbose, TString histoFileName)
: FairTask(name, iVerbose)
@@ -275,10 +280,20 @@ InitStatus CbmBbaAlignmentTask::Init()
}
for (int i = 0; i < fNtrackingStations + 1; i++) {
+ // set line colors
hResidualsAfterAlignmentX[i]->SetLineColor(kRed);
hResidualsAfterAlignmentY[i]->SetLineColor(kRed);
hPullsAfterAlignmentX[i]->SetLineColor(kRed);
hPullsAfterAlignmentY[i]->SetLineColor(kRed);
+ // set histograms to dynamic binning
+ hResidualsBeforeAlignmentX[i]->SetCanExtend(TH1::kXaxis);
+ hResidualsBeforeAlignmentY[i]->SetCanExtend(TH1::kXaxis);
+ hResidualsAfterAlignmentX[i]->SetCanExtend(TH1::kXaxis);
+ hResidualsAfterAlignmentY[i]->SetCanExtend(TH1::kXaxis);
+ hPullsBeforeAlignmentX[i]->SetCanExtend(TH1::kXaxis);
+ hPullsBeforeAlignmentY[i]->SetCanExtend(TH1::kXaxis);
+ hPullsAfterAlignmentX[i]->SetCanExtend(TH1::kXaxis);
+ hPullsAfterAlignmentY[i]->SetCanExtend(TH1::kXaxis);
}
gDirectory = curDirectory;
@@ -305,6 +320,8 @@ void CbmBbaAlignmentTask::Exec(Option_t* /*opt*/)
// select tracks for alignment and store them
unsigned num_rejected_fit = 0;
if (fTrackingMode == kSts && fInputStsTracks) {
+ int numPosMomentum = 0;
+ int numNegMomentum = 0;
fFitter.SetDefaultMomentumForMs(0.1);
@@ -337,10 +354,20 @@ void CbmBbaAlignmentTask::Exec(Option_t* /*opt*/)
if (t.fNstsHits < l1Par.GetNstationsActive(ca::EDetectorID::kSts)) continue;
if (!std::isfinite(par.GetQp())) continue;
if (fabs(par.GetQp()) > 1.) continue; // select tracks with min 1 GeV momentum
-
+ if (par.GetQp() > 0.) {
+ // if (numPosMomentum > numNegMomentum) {
+ // continue;
+ // }
+ numPosMomentum++;
+ }
+ else {
+ numNegMomentum++;
+ }
t.fAlignedTrack = t.fUnalignedTrack;
fTracks.push_back(t);
}
+ LOG(info) << "BBA: selected " << fTracks.size() << " tracks for alignment, " << numPosMomentum << " positive and "
+ << numNegMomentum << " negative momentum tracks";
}
else if (fTrackingMode == kMcbm && fInputGlobalTracks) {
@@ -484,6 +511,7 @@ void CbmBbaAlignmentTask::Exec(Option_t* /*opt*/)
}
+// This function destroys alignment if stations are defined as alignment bodies, do not use
void CbmBbaAlignmentTask::ConstrainStation(std::vector<double>& par, int iSta, int ixyz)
{
// set overall shifts of the station to zero
@@ -514,14 +542,14 @@ void CbmBbaAlignmentTask::ApplyConstraints(std::vector<double>& par)
ConstrainStation(par, fNtrackingStations - 1, 2);
}
-
+// This function destroys alignment if stations are defined as alignment bodies, do not use
void CbmBbaAlignmentTask::ApplyAlignment(const std::vector<double>& par)
{
// apply alignment parameters to the hits
std::vector<double> parConstrained = par;
- ApplyConstraints(parConstrained);
+ //ApplyConstraints(parConstrained);
for (TrackContainer& t : fTracks) {
@@ -784,26 +812,21 @@ void CbmBbaAlignmentTask::Finish()
bba::Parameter& p = par[ip];
p.SetActive(0);
p.SetValue(0.);
- p.SetBoundaryMin(-20.); // +-20 cm alignment range
- p.SetBoundaryMax(20.);
+ p.SetBoundaryMin(-2.);
+ p.SetBoundaryMax(2.);
p.SetStepMin(1.e-4);
p.SetStepMax(0.1);
p.SetStep(1.e-2);
}
// set active parameters
-
+ // fix first and last station
for (unsigned int ib = 0; ib < fAlignmentBodies.size(); ib++) {
auto& b = fAlignmentBodies[ib];
- if (b.fTrackingStation == 0) { // the first station
- //continue;
- }
- if (b.fTrackingStation == fNtrackingStations - 1) { // the last station
- //continue;
- }
- if (fTrackingMode == kSts && b.fTrackingStation == fNtrackingStations / 2) { // station in the middle for STS mode
+ if ((fTrackingMode == kSts && b.fTrackingStation == 0)
+ || (fTrackingMode == kSts && b.fTrackingStation == fNtrackingStations - 1)) {
par[3 * ib + 0].SetActive(0);
- par[3 * ib + 1].SetActive(1);
+ par[3 * ib + 1].SetActive(0);
par[3 * ib + 2].SetActive(0);
}
else {
@@ -811,8 +834,13 @@ void CbmBbaAlignmentTask::Finish()
par[3 * ib + 1].SetActive(1);
par[3 * ib + 2].SetActive(1);
}
+ // fix x in the middle station
+ if (fTrackingMode == kSts && b.fTrackingStation == fNtrackingStations / 2) {
+ par[3 * ib + 0].SetActive(0);
+ }
}
+
// set parameter ranges for different subsystems
for (const auto& s : fSensors) {
int ib = s.fAlignmentBody;
@@ -831,6 +859,8 @@ void CbmBbaAlignmentTask::Finish()
}
gRandom->SetSeed(1);
+ LOG(info) << "Simulated misalignment range: " << fSimulatedMisalignmentRange;
+ double d = fSimulatedMisalignmentRange;
for (int is = 0; is < fNalignmentBodies; is++) {
bba::Parameter& px = par[3 * is + 0];
@@ -838,7 +868,6 @@ void CbmBbaAlignmentTask::Finish()
bba::Parameter& pz = par[3 * is + 2];
// +- d cm misalignment
- double d = fSimulatedMisalignmentRange;
if (px.IsActive()) {
px.SetValue(gRandom->Uniform(2. * d) - d);
}
@@ -986,22 +1015,22 @@ void CbmBbaAlignmentTask::Finish()
continue;
}
const auto* n = el->GetPnode();
- LOG(info) << "Node: ";
- LOG(info) << n->GetName();
+ LOG(debug) << "Node: ";
+ LOG(debug) << n->GetName();
alignmentMatrices.insert(
AlignNode(n->GetName(), result[i * 3 + 0], result[i * 3 + 1], result[i * 3 + 2], 0., 0., 0.));
}
// for Trd stations
else if (s.fSystemId == ECbmModuleId::kTrd) {
- LOG(info) << "Node: ";
- LOG(info) << s.fNodePath;
+ LOG(debug) << "Node: ";
+ LOG(debug) << s.fNodePath;
alignmentMatrices.insert(
AlignNode(s.fNodePath, result[i * 3 + 0], result[i * 3 + 1], result[i * 3 + 2], 0., 0., 0.));
}
// for Trd2D station
else if (s.fSystemId == ECbmModuleId::kTrd2d) {
- LOG(info) << "Node: ";
- LOG(info) << s.fNodePath;
+ LOG(debug) << "Node: ";
+ LOG(debug) << s.fNodePath;
alignmentMatrices.insert(
AlignNode(s.fNodePath, result[i * 3 + 0], result[i * 3 + 1], result[i * 3 + 2], 0., 0., 0.));
}
@@ -1023,11 +1052,6 @@ void CbmBbaAlignmentTask::Finish()
misalignmentMatrixRootfile->Close();
}
-
- LOG(info) << " Cost function for the true parameters: " << fCostIdeal;
- LOG(info) << " Cost function for the initial parameters: " << fCostInitial;
- LOG(info) << " Cost function for the aligned parameters: " << costResult;
-
LOG(info) << "\n Misaligned parameters: ";
for (int is = 0; is < fNalignmentBodies; is++) {
const std::vector<double>& r = parMisaligned;
@@ -1044,16 +1068,10 @@ void CbmBbaAlignmentTask::Finish()
}
}
- LOG(info) << "\n";
-
- LOG(info) << " Cost function for the true parameters: " << fCostIdeal;
- LOG(info) << " Cost function for the initial parameters: " << fCostInitial;
- LOG(info) << " Cost function for the aligned parameters: " << costResult;
-
LOG(info) << "\n Alignment results, constrained: ";
{
std::vector<double> r = alignment.getResult();
- ApplyConstraints(r);
+ //ApplyConstraints(r);
for (int is = 0; is < fNalignmentBodies; is++) {
auto& b = fAlignmentBodies[is];
LOG(info) << "Body " << is << " sta " << b.fTrackingStation << ": x " << r[3 * is + 0] << " y " << r[3 * is + 1]
@@ -1128,6 +1146,7 @@ void CbmBbaAlignmentTask::Finish()
gDirectory = curr;
}
+
void CbmBbaAlignmentTask::WriteHistosCurFile(TObject* obj)
{
if (!obj->IsFolder())
diff --git a/reco/alignment/CbmBbaAlignmentTask.h b/reco/alignment/CbmBbaAlignmentTask.h
index f7bb9ffd35ad892df0b3c3b12c0d08bddd8186d4..0aa2f513198dec1a00790cfcefb61d62f9ee7c27 100644
--- a/reco/alignment/CbmBbaAlignmentTask.h
+++ b/reco/alignment/CbmBbaAlignmentTask.h
@@ -30,12 +30,7 @@ class TH1F;
///
/// an example of alignment using BBA package
-///
-/// you need to switch to the double precision in /algo/ca/CaSimdVc.h
-/// by uncommenting this line there:
-///
-/// typedef Vc::double_v fvec;
-///
+
class CbmBbaAlignmentTask : public FairTask {
public:
// Constructors/Destructors ---------
@@ -53,6 +48,7 @@ class CbmBbaAlignmentTask : public FairTask {
void SetStsTrackingMode() { fTrackingMode = TrackingMode::kSts; }
void SetSimulatedMisalignmentRange(double range) { fSimulatedMisalignmentRange = range; }
+ void SetRandomSeed(int seed) { fRandomSeed = seed; }
public:
enum TrackingMode
@@ -154,6 +150,8 @@ class CbmBbaAlignmentTask : public FairTask {
double fSimulatedMisalignmentRange{0.}; // misalignment range for simulated misalignment
+ int fRandomSeed{1};
+
double fChi2Total{0.};
long fNdfTotal{0};
long fFixedNdf{-1};