Bba: alignment task for the mcbm setup

macro/alignment/run_BbaAlignment_mcbm.C

+/* Copyright (C) 2020 Facility for Antiproton and Ion Research in Europe, Darmstadt
+   SPDX-License-Identifier: GPL-3.0-only
+   Authors: Dominik Smith [committer] */
+
+// --------------------------------------------------------------------------
+//
+// Macro for simulation & reconstruction QA
+//
+// The following naming conventions are assumed:
+// Raw data file:  [dataset].event.raw.root
+// Transport file: [dataset].tra.root
+// Parameter file:  [dataset].par.root
+// Reconstruction file: [dataset].rec.root
+//
+// S. Gorbunov 28/09/2020
+//
+// --------------------------------------------------------------------------
+
+// Includes needed for IDE
+#if !defined(__CLING__)
+#include "CbmDefs.h"
+#include "CbmMCDataManager.h"
+#include "CbmMuchTransportQa.h"
+#include "CbmQaManager.h"
+#include "CbmSetup.h"
+
+#include <FairFileSource.h>
+#include <FairMonitor.h>
+#include <FairParAsciiFileIo.h>
+#include <FairParRootFileIo.h>
+#include <FairRootFileSink.h>
+#include <FairRunAna.h>
+#include <FairRuntimeDb.h>
+#include <FairSystemInfo.h>
+
+#include <TStopwatch.h>
+#endif
+
+void run_BbaAlignment_mcbm(Int_t nEvents = 0, TString dataset = "data/mcbm_beam_2020_03_test",
+                           TString setupName = "mcbm_beam_2020_03", Bool_t bUseMC = kTRUE, const TString& config = "")
+{
+
+  // ========================================================================
+  //          Adjust this part according to your requirements
+
+  // -----   Logger settings   ----------------------------------------------
+  FairLogger::GetLogger()->SetLogScreenLevel("INFO");
+  FairLogger::GetLogger()->SetLogVerbosityLevel("LOW");
+  // ------------------------------------------------------------------------
+
+  // -----   Environment   --------------------------------------------------
+  int verbose      = 6;                              // verbose level
+  TString myName   = "mcbm_qa";                      // this macro's name for screen output
+  TString srcDir   = gSystem->Getenv("VMCWORKDIR");  // top source directory
+  TString qaConfig = (config.Length() ? config : srcDir + "/macro/qa/configs/qa_tasks_config_mcbm.yaml");
+  // NOTE: SZh 28.07.2024: config can depend from the setup
+  // ------------------------------------------------------------------------
+
+  // -----   In- and output file names   ------------------------------------
+  TString rawFile  = dataset + ".event.raw.root";
+  TString traFile  = dataset + ".tra.root";
+  TString parFile  = dataset + ".par.root";
+  TString recFile  = dataset + ".rec.root";
+  TString sinkFile = dataset + ".qa.root";
+  // ------------------------------------------------------------------------
+
+  // -----   Load the geometry setup   -------------------------------------
+  std::cout << '\n';
+  TString setupFile  = srcDir + "/geometry/setup/setup_" + setupName + ".C";
+  TString setupFunct = "setup_";
+  setupFunct         = setupFunct + setupName + "()";
+  std::cout << "-I- " << myName << ": Loading macro " << setupFile << '\n';
+  gROOT->LoadMacro(setupFile);
+  gROOT->ProcessLine(setupFunct);
+  CbmSetup* setup = CbmSetup::Instance();
+  //  setup->RemoveModule(ECbmModuleId::kTrd);
+  // ------------------------------------------------------------------------
+
+  // -----   Some global switches   -----------------------------------------
+  //if (setupName == "mcbm_beam_2022_05_23_nickel") { setup->RemoveModule(ECbmModuleId::kMuch); }
+
+  //bool eventBased = !sEvBuildRaw.IsNull();
+  bool bUseMvd  = setup->IsActive(ECbmModuleId::kMvd);
+  bool bUseSts  = setup->IsActive(ECbmModuleId::kSts);
+  bool bUseRich = setup->IsActive(ECbmModuleId::kRich);
+  bool bUseMuch = setup->IsActive(ECbmModuleId::kMuch);
+  bool bUseTrd  = setup->IsActive(ECbmModuleId::kTrd);
+  bool bUseTof  = setup->IsActive(ECbmModuleId::kTof);
+  bool bUsePsd  = setup->IsActive(ECbmModuleId::kPsd);
+  std::cout << "  MVD: " << (bUseMvd ? "ON" : "OFF") << '\n';
+  std::cout << "  STS: " << (bUseSts ? "ON" : "OFF") << '\n';
+  std::cout << "  RICH: " << (bUseRich ? "ON" : "OFF") << '\n';
+  std::cout << "  MUCH: " << (bUseMuch ? "ON" : "OFF") << '\n';
+  std::cout << "  TRD: " << (bUseTrd ? "ON" : "OFF") << '\n';
+  std::cout << "  TOF: " << (bUseTof ? "ON" : "OFF") << '\n';
+  std::cout << "  PSD: " << (bUsePsd ? "ON" : "OFF") << '\n';
+  // ------------------------------------------------------------------------
+
+  // -----   Parameter files as input to the runtime database   -------------
+  std::cout << '\n';
+  std::cout << "-I- " << myName << ": Defining paramete files\n";
+  TList* parFileList = new TList();
+  TString geoTag;
+
+  // - MUCH digitisation parameters
+  TString muchParFile {};
+  if (setup->GetGeoTag(ECbmModuleId::kMuch, geoTag)) {
+    bool mcbmFlag = geoTag.Contains("mcbm", TString::kIgnoreCase);
+    muchParFile   = srcDir + "/parameters/much/much_";
+    muchParFile += (mcbmFlag) ? geoTag : geoTag(0, 4);
+    muchParFile += "_digi_sector.root";
+    {  // init geometry from the file
+      TFile* f            = new TFile(muchParFile, "R");
+      TObjArray* stations = f->Get<TObjArray>("stations");
+      assert(stations);
+      CbmMuchGeoScheme::Instance()->Init(stations, mcbmFlag);
+    }
+  }
+
+  // - TRD digitisation parameters
+  if (setup->GetGeoTag(ECbmModuleId::kTrd, geoTag)) {
+    const Char_t* npar[4] = {"asic", "digi", "gas", "gain"};
+    TObjString* trdParFile(NULL);
+    for (Int_t i(0); i < 4; i++) {
+      trdParFile = new TObjString(srcDir + "/parameters/trd/trd_" + geoTag + "." + npar[i] + ".par");
+      parFileList->Add(trdParFile);
+      std::cout << "-I- " << myName << ": Using parameter file " << trdParFile->GetString() << std::endl;
+    }
+  }
+
+  // - TOF digitisation parameters
+  if (setup->GetGeoTag(ECbmModuleId::kTof, geoTag)) {
+    TObjString* tofBdfFile = new TObjString(srcDir + "/parameters/tof/tof_" + geoTag + ".digibdf.par");
+    parFileList->Add(tofBdfFile);
+    std::cout << "-I- " << myName << ": Using parameter file " << tofBdfFile->GetString() << std::endl;
+  }
+  // ------------------------------------------------------------------------
+
+  // In general, the following parts need not be touched
+  // ========================================================================
+
+  // -----   Timer   --------------------------------------------------------
+  TStopwatch timer;
+  timer.Start();
+  // ------------------------------------------------------------------------
+
+  // ----    Debug option   -------------------------------------------------
+  gDebug = 0;
+  // ------------------------------------------------------------------------
+
+  // -----   FairRunAna   ---------------------------------------------------
+  FairFileSource* inputSource = new FairFileSource(rawFile);
+  inputSource->AddFriend(traFile);
+  inputSource->AddFriend(recFile);
+
+  FairRunAna* run = new FairRunAna();
+  run->SetSource(inputSource);
+  run->SetGenerateRunInfo(kFALSE);
+
+
+  FairRootFileSink* sink = new FairRootFileSink(sinkFile);
+  run->SetSink(sink);
+
+  TString monitorFile {sinkFile};
+  monitorFile.ReplaceAll("qa", "qa.monitor");
+  FairMonitor::GetMonitor()->EnableMonitor(kTRUE, monitorFile);
+  // ------------------------------------------------------------------------
+
+  // -----   MCDataManager (legacy mode)  -----------------------------------
+  if (bUseMC) {
+    std::cout << "-I- " << myName << ": Adding MC manager and MC to reco matching tasks\n";
+
+    auto* mcManager = new CbmMCDataManager("MCDataManager", 1);
+    mcManager->AddFile(traFile);
+    run->AddTask(mcManager);
+
+    auto* matchRecoToMC = new CbmMatchRecoToMC();
+    // NOTE: Matching is suppressed, if there are hit and cluster matches already in the tree. If there
+    //       are no hit matches, they are produced on this stage.
+    matchRecoToMC->SuppressHitReMatching();
+    run->AddTask(matchRecoToMC);
+  }
+  // ------------------------------------------------------------------------
+
+  // ----- Tracking detector interface --------------------------------------
+  run->AddTask(new CbmTrackingDetectorInterfaceInit());
+  // ------------------------------------------------------------------------
+
+
+  // ----- CA tracking QA ---------------------------------------------------
+  // Kalman Filter (currently needed to access the magnetic filed, to be
+  // removed soon)
+  run->AddTask(new CbmKF());
+
+  // ------------------------------------------------------------------------
+
+  // TODO: read tracking parameters from the file like CaOutputQa does
+  //TODO:  instead of initializing the tracker
+  //TString caParFile = recFile;
+  //caParFile.ReplaceAll(".root", ".L1Parameters.dat");
+  //auto* pCaOutputQa = new cbm::ca::OutputQa(verbose, bUseMC);
+  //pCaOutputQa->ReadParameters(caParFile.Data());
+
+  // L1 CA track finder setup
+  auto l1 = new CbmL1("CA");
+  l1->SetMcbmMode();
+  l1->DisableTrackingStation(cbm::algo::ca::EDetectorID::kMuch, 0);
+  l1->DisableTrackingStation(cbm::algo::ca::EDetectorID::kMuch, 1);
+  l1->DisableTrackingStation(cbm::algo::ca::EDetectorID::kMuch, 2);
+
+  // User configuration example for CA:
+  //l1->SetConfigUser(srcDir + "/macro/L1/configs/ca_params_user_example.yaml");
+  run->AddTask(l1);
+
+  // ----- BBA alignment   --------------------------------------------
+  CbmBbaAlignmentTask* alignemnt = new CbmBbaAlignmentTask();
+  //CbmKFTrackQa* alignemnt = new CbmKFTrackQa();
+  run->AddTask(alignemnt);
+
+
+  // -----  Parameter database   --------------------------------------------
+  std::cout << std::endl << std::endl;
+  std::cout << "-I- " << myName << ": Set runtime DB" << std::endl;
+  FairRuntimeDb* rtdb       = run->GetRuntimeDb();
+  FairParRootFileIo* parIo1 = new FairParRootFileIo();
+  parIo1->open(parFile.Data(), "in");
+  rtdb->setFirstInput(parIo1);
+  if (!parFileList->IsEmpty()) {
+    FairParAsciiFileIo* parIo2 = new FairParAsciiFileIo();
+    parIo2->open(parFileList, "in");
+    rtdb->setSecondInput(parIo2);
+  }
+  rtdb->print();
+  // ------------------------------------------------------------------------
+
+  // -----   Run initialisation   -------------------------------------------
+  std::cout << std::endl;
+  std::cout << "-I- " << myName << ": Initialise run" << std::endl;
+  run->Init();
+
+  // -----   Start run   ----------------------------------------------------
+  std::cout << std::endl << std::endl;
+  std::cout << "-I- " << myName << ": Starting run" << std::endl;
+  run->Run(0, nEvents);
+  // ------------------------------------------------------------------------
+
+  // -----   Finish   -------------------------------------------------------
+  timer.Stop();
+
+  FairMonitor::GetMonitor()->Print();
+  Double_t rtime = timer.RealTime();
+  Double_t ctime = timer.CpuTime();
+  std::cout << std::endl << std::endl;
+  std::cout << "Macro finished successfully." << std::endl;
+  std::cout << "Output file is " << sinkFile << std::endl;
+  std::cout << "Parameter file is " << parFile << std::endl;
+  std::cout << "Real time " << rtime << " s, CPU time " << ctime << " s" << std::endl;
+  std::cout << std::endl;
+  std::cout << " Test passed" << std::endl;
+  std::cout << " All ok " << std::endl;
+  // ------------------------------------------------------------------------
+
+  // -----   Resource monitoring   ------------------------------------------
+  // Extract the maximal used memory an add is as Dart measurement
+  // This line is filtered by CTest and the value send to CDash
+  FairSystemInfo sysInfo;
+  Float_t maxMemory = sysInfo.GetMaxMemory();
+  std::cout << "<DartMeasurement name=\"MaxMemory\" type=\"numeric/double\">";
+  std::cout << maxMemory;
+  std::cout << "</DartMeasurement>" << std::endl;
+
+  Float_t cpuUsage = ctime / rtime;
+  std::cout << "<DartMeasurement name=\"CpuLoad\" type=\"numeric/double\">";
+  std::cout << cpuUsage;
+  std::cout << "</DartMeasurement>" << std::endl;
+  // ------------------------------------------------------------------------
+
+  // -----   Function needed for CTest runtime dependency   -----------------
+  RemoveGeoManager();
+  // ------------------------------------------------------------------------
+}

reco/alignment/CMakeLists.txt

@@ -17,6 +17,12 @@ set(PUBLIC_DEPENDENCIES
 )
 
 set(PRIVATE_DEPENDENCIES
+  CbmMvdBase
+  CbmStsBase
+  CbmMuchBase
+  CbmTrdBase
+  CbmTofBase
+  L1
   KF
   bba::library
   CbmStsBase

reco/alignment/CbmBbaAlignmentTask.cxx

@@ -12,7 +12,6 @@
 // Cbm Headers ----------------------
 #include "CbmBbaAlignmentTask.h"
 
-#include "CbmL1PFFitter.h"
 
 // ROOT headers
 
@@ -24,11 +23,21 @@
 #include "TRandom.h"
 
 //
-
+#include "CbmCaTimeSliceReader.h"
+#include "CbmGlobalTrack.h"
+#include "CbmL1.h"
+#include "CbmMuchTrack.h"
+#include "CbmMuchTrackingInterface.h"
 #include "CbmMvdHit.h"
+#include "CbmMvdTrackingInterface.h"
 #include "CbmStsHit.h"
 #include "CbmStsSetup.h"
 #include "CbmStsTrack.h"
+#include "CbmStsTrackingInterface.h"
+#include "CbmTofTrack.h"
+#include "CbmTofTrackingInterface.h"
+#include "CbmTrdTrack.h"
+#include "CbmTrdTrackingInterface.h"
 
 #include "bba/BBA.h"
 
@@ -36,6 +45,30 @@
 
 #include <iostream>
 
+namespace
+{
+  using namespace cbm::algo;
+}
+
+void CbmBbaAlignmentTask::TrackContainer::MakeConsistent()
+{
+  fNmvdHits  = 0;
+  fNstsHits  = 0;
+  fNmuchHits = 0;
+  fNtrdHits  = 0;
+  fNtofHits  = 0;
+  for (const auto& n : fUnalignedTrack.fNodes) {
+    switch (n.fSystemId) {
+      case ECbmModuleId::kMvd: fNmvdHits++; break;
+      case ECbmModuleId::kSts: fNstsHits++; break;
+      case ECbmModuleId::kMuch: fNmuchHits++; break;
+      case ECbmModuleId::kTrd: fNtrdHits++; break;
+      case ECbmModuleId::kTof: fNtofHits++; break;
+      default: break;
+    }
+  }
+}
+
 
 CbmBbaAlignmentTask::CbmBbaAlignmentTask(const char* name, Int_t iVerbose, TString histoFileName)
   : FairTask(name, iVerbose)
@@ -66,6 +99,9 @@ CbmBbaAlignmentTask::~CbmBbaAlignmentTask() {}
 
 InitStatus CbmBbaAlignmentTask::Init()
 {
+
+  fFitter.Init();
+
   //Get ROOT Manager
   FairRootManager* ioman = FairRootManager::Instance();
 
@@ -74,35 +110,45 @@ InitStatus CbmBbaAlignmentTask::Init()
     return kERROR;
   }
 
-  // Get hits
+  // Get global tracks & matches
+
+  if (TrackingMode::kSts != fTrackingMode) {
+
+    fInputGlobalTracks = static_cast<TClonesArray*>(ioman->GetObject("GlobalTrack"));
+    if (!fInputGlobalTracks) {
+      LOG(error) << "CbmBbaAlignmentTask::Init: Global track array not found!";
+      return kERROR;
+    }
+
+    fInputGlobalTrackMatches = static_cast<TClonesArray*>(ioman->GetObject("GlobalTrackMatch"));
 
-  fInputMvdHits = static_cast<TClonesArray*>(ioman->GetObject("MvdHit"));
-  fInputStsHits = static_cast<TClonesArray*>(ioman->GetObject("StsHit"));
+    if (!fInputGlobalTrackMatches) {
+      LOG(error) << "CbmBbaAlignmentTask::Init: Global track matches array not found!";
+      //return kERROR;
+    }
+  }
 
-  // Get sts tracks
   fInputStsTracks = static_cast<TClonesArray*>(ioman->GetObject("StsTrack"));
   if (!fInputStsTracks) {
-    LOG(error) << "CbmBbaAlignmentTask::Init: Sts track-array not found!";
+    LOG(error) << "CbmBbaAlignmentTask::Init: Sts track array not found!";
     return kERROR;
   }
 
-  // MC track match
+  fInputStsTrackMatches = static_cast<TClonesArray*>(ioman->GetObject("StsTrackMatch"));
+
+  if (!fInputStsTrackMatches) {
+    LOG(error) << "CbmBbaAlignmentTask::Init: Sts track matches array not found!";
+    //return kERROR;
+  }
+
+  // MC tracks
   fInputMcTracks = static_cast<TClonesArray*>(ioman->GetObject("MCTrack"));
   if (!fInputMcTracks) {
     Warning("CbmBbaAlignmentTask::Init", "mc track array not found!");
     return kERROR;
   }
 
-  // Track match
-  fInputStsTrackMatches = static_cast<TClonesArray*>(ioman->GetObject("StsTrackMatch"));
-  if (fInputStsTrackMatches == 0) {
-    LOG(error) << "CbmBbaAlignmentTask::Init: track match array not found!";
-    return kERROR;
-  }
-
   fTracks.clear();
-  fMvdHits.clear();
-  fStsHits.clear();
 
   return kSUCCESS;
 }
@@ -116,106 +162,110 @@ void CbmBbaAlignmentTask::Exec(Option_t* /*opt*/)
 
   if (static_cast<int>(fTracks.size()) >= fMaxNtracks) { return; }
 
-  // select STS tracks for alignment and store them
-
-  for (int iTr = 0; iTr < fInputStsTracks->GetEntriesFast(); iTr++) {
-
-    if (static_cast<int>(fTracks.size()) >= fMaxNtracks) { break; }
-
-    CbmStsTrack* stsTrack = static_cast<CbmStsTrack*>(fInputStsTracks->At(iTr));
-
-    if (stsTrack->GetNofStsHits() < 8) continue;
-    const auto* par = stsTrack->GetParamFirst();
-    if (!std::isfinite(par->GetQp())) continue;
-    if (fabs(par->GetQp()) > 1.) continue;  // select tracks with min 1 GeV momentum
-
-    CbmStsTrack track;
-
-    // copy track parameters
-
-    track.SetParamFirst(stsTrack->GetParamFirst());
-    track.SetParamLast(stsTrack->GetParamLast());
-    track.SetChiSq(stsTrack->GetChiSq());
-    track.SetNDF(stsTrack->GetNDF());
-
-    // copy hits to the local arrays and set the track hit indices correspondingly
-
-    for (int ih = 0; ih < stsTrack->GetNofMvdHits(); ih++) {
-      assert(fInputMvdHits);
-      int hitIndex  = stsTrack->GetMvdHitIndex(ih);
-      CbmMvdHit hit = *dynamic_cast<CbmMvdHit*>(fInputMvdHits->At(hitIndex));
-      track.AddMvdHit(fMvdHits.size());
-      fMvdHits.push_back(hit);
+  ca::Framework* l1           = CbmL1::Instance()->fpAlgo;
+  const ca::Parameters& l1Par = l1->GetParameters();
+
+  static int statGlobalTracks = 0;
+  statGlobalTracks += fInputGlobalTracks->GetEntriesFast();
+  static int statGlobalTracks1 = 0;
+  static int statGlobalTracks2 = 0;
+  // select tracks for alignment and store them
+  if (fTrackingMode == kSts && fInputStsTracks) {
+    for (int iTr = 0; iTr < fInputStsTracks->GetEntriesFast(); iTr++) {
+      if (static_cast<int>(fTracks.size()) >= fMaxNtracks) { break; }
+      const CbmStsTrack* stsTrack = dynamic_cast<CbmStsTrack*>(fInputStsTracks->At(iTr));
+      if (!stsTrack) continue;
+      TrackContainer t;
+      if (!fFitter.CreateMvdStsTrack(t.fUnalignedTrack, *stsTrack)) continue;
+      t.MakeConsistent();
+      fFitter.FitTrack(t.fUnalignedTrack);
+      t.fAlignedTrack = t.fUnalignedTrack;
+      if (t.fNstsHits < l1Par.GetNstationsActive(ca::EDetectorID::kSts)) continue;
+      const auto& par = t.fUnalignedTrack.fNodes[t.fUnalignedTrack.fFirstHitNode].fTrack;
+      if (!std::isfinite(par.GetQp()[0])) continue;
+      if (fabs(par.GetQp()[0]) > 1.) continue;  // select tracks with min 1 GeV momentum
+      fTracks.push_back(t);
     }
-
-    for (int ih = 0; ih < stsTrack->GetNofStsHits(); ih++) {
-      assert(fInputStsHits);
-      int hitIndex  = stsTrack->GetStsHitIndex(ih);
-      CbmStsHit hit = *dynamic_cast<CbmStsHit*>(fInputStsHits->At(hitIndex));
-      track.AddStsHit(fStsHits.size());
-      fStsHits.push_back(hit);
+  }
+  else if (fTrackingMode == kMcbm && fInputGlobalTracks) {
+    for (int iTr = 0; iTr < fInputGlobalTracks->GetEntriesFast(); iTr++) {
+      if (static_cast<int>(fTracks.size()) >= fMaxNtracks) { break; }
+      const CbmGlobalTrack* globalTrack = dynamic_cast<const CbmGlobalTrack*>(fInputGlobalTracks->At(iTr));
+      if (!globalTrack) continue;
+      statGlobalTracks1++;
+      TrackContainer t;
+      if (!fFitter.CreateGlobalTrack(t.fUnalignedTrack, *globalTrack)) {
+        LOG(error) << "BBA: can not create a global track for the fit! ";
+        exit(0);
+        continue;
+      }
+      statGlobalTracks2++;
+      t.MakeConsistent();
+      fFitter.FitTrack(t.fUnalignedTrack);
+      t.fAlignedTrack = t.fUnalignedTrack;
+      //if (t.fNstsHits < 1) continue;
+      //if (t.fNtrdHits < 2) continue;
+      fTracks.push_back(t);
     }
+  }
 
-    fTracks.push_back(track);
-
-  }  // tracks
+  LOG(info) << "Bba: " << fInputGlobalTracks->GetEntriesFast() << " global tracks in event, " << statGlobalTracks
+            << " global tracks in total, " << fTracks.size() << " tracks selected for alignment";
 }
 
+
 void CbmBbaAlignmentTask::ApplyAlignment(const std::vector<double>& par)
 {
   // apply alignment parameters to the hits
-  for (unsigned int ih = 0; ih < fStsHits.size(); ih++) {
-
-    const CbmStsHit& hitOrig = fStsHits[ih];
 
-    int iStation = CbmStsSetup::Instance()->GetStationNumber(hitOrig.GetAddress());
+  for (TrackContainer& t : fTracks) {
 
-    assert(iStation >= 0 && iStation < CbmStsSetup::Instance()->GetNofStations());
+    for (unsigned int in = 0; in < t.fUnalignedTrack.fNodes.size(); in++) {
+      CbmKFTrackFitter::FitNode& node        = t.fUnalignedTrack.fNodes[in];
+      CbmKFTrackFitter::FitNode& nodeAligned = t.fAlignedTrack.fNodes[in];
+      int iAlignmentModule                   = node.fReference1;
+      if (iAlignmentModule < 0) continue;
+      assert(nodeAligned.fReference1 == iAlignmentModule);
 
-    double dx = par[3 * iStation + 0];
-    double dy = par[3 * iStation + 1];
-    double dz = par[3 * iStation + 2];
+      double dx = par[3 * iAlignmentModule + 0];
+      double dy = par[3 * iAlignmentModule + 1];
+      double dz = par[3 * iAlignmentModule + 2];
 
-    CbmStsHit& hitAligned = fStsHitsAligned[ih];
-
-    hitAligned.SetX(hitOrig.GetX() + dx);
-    hitAligned.SetY(hitOrig.GetY() + dy);
-    hitAligned.SetZ(hitOrig.GetZ() + dz);
+      nodeAligned.fMxy.SetX(node.fMxy.X() + ca::fvec(dx));
+      nodeAligned.fMxy.SetY(node.fMxy.Y() + ca::fvec(dy));
+      nodeAligned.fZ = node.fZ + dz;
+    }
   }
 }
 
+
 double CbmBbaAlignmentTask::CostFunction(const std::vector<double>& par)
 {
-
   // apply new parameters to the hits
 
   ApplyAlignment(par);
 
-  auto newTracks = fTracks;
-
-  CbmL1PFFitter fitter;
-  fitter.Fit(newTracks, fMvdHits, fStsHitsAligned, fTrackPdg);
-
-  double chi2Total  = 0;
-  long int ndfTotal = 1;
+  fChi2Total = 0.;
+  fNdfTotal  = 0;
 
-  int nGoodTracks = 0;
-  for (unsigned int iTr = 0; iTr < newTracks.size(); iTr++) {
-    if (!std::isfinite(newTracks[iTr].GetChiSq())) continue;
-    if (newTracks[iTr].GetChiSq() < 0.) continue;
-    if (newTracks[iTr].GetNDF() < 0.) continue;
-
-    nGoodTracks++;
-    chi2Total += newTracks[iTr].GetChiSq();
-    ndfTotal += newTracks[iTr].GetNDF();
+  for (auto& t : fTracks) {
+    fFitter.FitTrack(t.fAlignedTrack);
+    double chi2 = t.fAlignedTrack.fNodes[t.fAlignedTrack.fFirstHitNode].fTrack.GetChiSq()[0];
+    double ndf  = t.fAlignedTrack.fNodes[t.fAlignedTrack.fFirstHitNode].fTrack.GetNdf()[0];
+    if (!std::isfinite(chi2) || chi2 <= 0. || ndf <= 0.) {
+      chi2 = 0.;
+      ndf  = 0.;
+    }
+    fChi2Total += chi2;
+    fNdfTotal += (int) ndf;
   }
 
-  double cost = chi2Total / ndfTotal;
-  LOG(info) << "BBA: cost function:  n tracks " << nGoodTracks << ", cost " << cost
+  double cost = fChi2Total / (fNdfTotal + 1);
+  LOG(info) << "BBA: cost function:  ndf " << fNdfTotal << ", cost " << cost
             << ", diff to ideal cost: " << cost - fCostIdeal;
+
+  if (fFixedNdf > 0 && fNdfTotal != fFixedNdf) { cost = 1.e30; }
   return cost;
-  if (nGoodTracks == static_cast<int>(fTracks.size())) { return cost; }
-  return 1.e30;
 }
 
 void CbmBbaAlignmentTask::Finish()
@@ -226,15 +276,22 @@ void CbmBbaAlignmentTask::Finish()
 
   LOG(info) << "BBA: start the alignment procedure with " << fTracks.size() << " tracks ...";
 
-  // init auxiliary arrays
+  ca::Framework* l1 = CbmL1::Instance()->fpAlgo;
 
-  fStsHitsAligned = fStsHits;
-  fTrackPdg.resize(fTracks.size(), 211);  // pion PDG
+  const ca::Parameters& l1Param = l1->GetParameters();
 
   // init alignmemt module
 
-  int nStsStations = CbmStsSetup::Instance()->GetNofStations();
-  int nParameters  = 3 * nStsStations;  // x,yz
+  int nStations = l1Param.GetNstationsActive();
+
+  for (auto& t : fTracks) {
+    for (auto& n : t.fUnalignedTrack.fNodes) {
+      n.fReference1 = n.fMaterialLayer;  // material layer is currently equal to the active tracking station
+    }
+    t.fAlignedTrack = t.fUnalignedTrack;
+  }
+
+  int nParameters = 3 * nStations;  // x, y, z
 
   std::vector<bba::Parameter> par(nParameters);
 
@@ -249,48 +306,25 @@ void CbmBbaAlignmentTask::Finish()
     p.SetStep(10.e-4);
   }
 
-  par[3 * (nStsStations - 1) + 0].SetActive(0);  // fix the last station
-  par[3 * (nStsStations - 1) + 1].SetActive(0);
-  par[3 * (nStsStations - 1) + 2].SetActive(0);
-
-  //par[3 * 1 + 1].SetActive(1);
-
-  par[3 * 0 + 0].SetActive(1);
-  par[3 * 0 + 1].SetActive(1);
-  par[3 * 0 + 2].SetActive(1);
-
-  par[3 * 1 + 0].SetActive(1);
-  par[3 * 1 + 1].SetActive(1);
-  par[3 * 1 + 2].SetActive(1);
-
-  par[3 * 2 + 0].SetActive(1);
-  par[3 * 2 + 1].SetActive(1);
-  par[3 * 2 + 2].SetActive(1);
-
-  par[3 * 3 + 0].SetActive(1);
-  par[3 * 3 + 1].SetActive(1);
-  par[3 * 3 + 2].SetActive(1);
+  par[0].SetActive(0);  // fix the first station
+  par[1].SetActive(0);
+  par[2].SetActive(0);
 
-  par[3 * 4 + 0].SetActive(1);
-  par[3 * 4 + 1].SetActive(1);
-  par[3 * 4 + 2].SetActive(1);
+  par[3 * (nStations - 1) + 0].SetActive(0);  // fix the last station
+  par[3 * (nStations - 1) + 1].SetActive(0);
+  par[3 * (nStations - 1) + 2].SetActive(0);
 
-  par[3 * 5 + 0].SetActive(1);
-  par[3 * 5 + 1].SetActive(1);
-  par[3 * 5 + 2].SetActive(1);
-
-  par[3 * 6 + 0].SetActive(1);
-  par[3 * 6 + 1].SetActive(1);
-  par[3 * 6 + 2].SetActive(1);
-
-  par[3 * 7 + 0].SetActive(0);
-  par[3 * 7 + 1].SetActive(0);
-  par[3 * 7 + 2].SetActive(0);
+  // set active parameters for other stations
 
+  for (int i = 1; i < nStations - 1; i++) {
+    par[3 * i + 0].SetActive(1);
+    par[3 * i + 1].SetActive(1);
+    par[3 * i + 2].SetActive(1);
+  }
 
   gRandom->SetSeed(1);
 
-  for (int is = 0; is < nStsStations; is++) {
+  for (int is = 0; is < nStations; is++) {
     bba::Parameter& px = par[3 * is + 0];
     bba::Parameter& py = par[3 * is + 1];
     bba::Parameter& pz = par[3 * is + 2];
@@ -327,33 +361,19 @@ void CbmBbaAlignmentTask::Finish()
   }
 
   {
-    ApplyAlignment(parMisaligned);
-
-    CbmL1PFFitter fitter;
-    auto oldTracks = fTracks;
-    auto newTracks = fTracks;
-    fitter.Fit(newTracks, fMvdHits, fStsHitsAligned, fTrackPdg);
-
-    double chi2Total  = 0;
-    long int ndfTotal = 1;
-
-    int nGoodTracks = 0;
+    fCostInitial   = CostFunction(parMisaligned);
+    auto tmpTracks = fTracks;
     fTracks.clear();
-    for (unsigned int iTr = 0; iTr < newTracks.size(); iTr++) {
-      if (!std::isfinite(newTracks[iTr].GetChiSq())) continue;
-      if (newTracks[iTr].GetChiSq() < 0.) continue;
-      if (newTracks[iTr].GetNDF() < 0.) continue;
-      fTracks.push_back(oldTracks[iTr]);
-      nGoodTracks++;
-      chi2Total += newTracks[iTr].GetChiSq();
-      ndfTotal += newTracks[iTr].GetNDF();
+    for (auto& t : tmpTracks) {
+      if (t.fAlignedTrack.fNodes[t.fAlignedTrack.fFirstHitNode].fTrack.GetNdf()[0] > 0) { fTracks.push_back(t); }
     }
-
-    LOG(info) << "Initial nTracks " << nGoodTracks << " chi2/ndf " << chi2Total / ndfTotal;
+    fFixedNdf = fNdfTotal;
   }
 
   fCostIdeal = CostFunction(parAligned);
 
+
+  LOG(info) << "Initial nTracks: " << fTracks.size() << " cost " << fCostInitial;
   LOG(info) << " cost function for the true parameters: " << fCostIdeal;
 
   alignment.align();
@@ -361,14 +381,14 @@ void CbmBbaAlignmentTask::Finish()
   LOG(info) << " cost function for the true parameters: " << fCostIdeal;
 
   LOG(info) << " Misaligned parameters: ";
-  for (int is = 0; is < nStsStations; is++) {
+  for (int is = 0; is < nStations; is++) {
     const std::vector<double>& r = parMisaligned;
     LOG(info) << "Sts station " << is << ": x " << r[3 * is + 0] << " y " << r[3 * is + 1] << " z " << r[3 * is + 2];
   }
 
   LOG(info) << " Alignment results: ";
 
-  for (int is = 0; is < nStsStations; is++) {
+  for (int is = 0; is < nStations; is++) {
     const std::vector<double>& r = alignment.getResult();
     LOG(info) << "Sts station " << is << ": x " << r[3 * is + 0] << " y " << r[3 * is + 1] << " z " << r[3 * is + 2];
   }

reco/alignment/CbmBbaAlignmentTask.h

@@ -12,13 +12,12 @@
 #ifndef CbmBbaAlignmentTask_H
 #define CbmBbaAlignmentTask_H
 
-#include "CbmMvdHit.h"
-#include "CbmStsHit.h"
-#include "CbmStsTrack.h"
+
+#include "CbmKFTrackFitter.h"
 
 #include "FairTask.h"
 
-#include <TString.h>
+#include "TString.h"
 
 #include <vector>
 
@@ -29,6 +28,7 @@ class TFile;
 class TDirectory;
 class TH1F;
 
+
 ///
 /// an example of alignment using BBA package
 ///
@@ -50,35 +50,57 @@ public:
   void Exec(Option_t* opt);
   void Finish();
 
+  void SetMcbmTrackingMode() { fTrackingMode = TrackingMode::kMcbm; }
+  void SetStsTrackingMode() { fTrackingMode = TrackingMode::kSts; }
+
+public:
+  enum TrackingMode
+  {
+    kSts,
+    kMcbm
+  };
+
+  /// information about a track
+  /// aligned and unaligned hits are stored in the corresponding CbmKFTrackFitter::Track objects
+  struct TrackContainer {
+
+    CbmKFTrackFitter::Track fUnalignedTrack;  // track before alignment
+    CbmKFTrackFitter::Track fAlignedTrack;    // track after alignment
+
+    int fNmvdHits {0};   // number of MVD hits
+    int fNstsHits {0};   // number of STS hits
+    int fNmuchHits {0};  // number of MUCH hits
+    int fNtrdHits {0};   // number of TRD hits
+    int fNtofHits {0};   // number of TOF hits
+
+    void MakeConsistent();
+  };
+
 private:
   const CbmBbaAlignmentTask& operator=(const CbmBbaAlignmentTask&);
   CbmBbaAlignmentTask(const CbmBbaAlignmentTask&);
 
   void WriteHistosCurFile(TObject* obj);
-  int GetHistoIndex(int pdg);
 
   void ApplyAlignment(const std::vector<double>& par);
 
   double CostFunction(const std::vector<double>& par);
 
-  //input data arrays
-  TClonesArray* fInputMvdHits {nullptr};
-  TClonesArray* fInputStsHits {nullptr};
-  TClonesArray* fInputStsTracks {nullptr};
+  TrackingMode fTrackingMode = TrackingMode::kMcbm;
 
-  TClonesArray* fInputMcTracks {nullptr};         // Mc info for debugging
-  TClonesArray* fInputStsTrackMatches {nullptr};  // Mc info for debugging
+  // input data arrays
 
-  // collection of selected tracks and hits
-  std::vector<CbmStsTrack> fTracks;
-  std::vector<CbmMvdHit> fMvdHits;
-  std::vector<CbmStsHit> fStsHits;
+  TClonesArray* fInputGlobalTracks {nullptr};
+  TClonesArray* fInputStsTracks {nullptr};
 
-  // temporary array with aligned hits
-  std::vector<CbmStsHit> fStsHitsAligned;
+  TClonesArray* fInputMcTracks {nullptr};            // Mc info for debugging
+  TClonesArray* fInputGlobalTrackMatches {nullptr};  // Mc info for debugging
+  TClonesArray* fInputStsTrackMatches {nullptr};     // Mc info for debugging
 
-  // array with pdg hypothesis for tracks
-  std::vector<int> fTrackPdg;
+  CbmKFTrackFitter fFitter;
+
+  // collection of selected tracks and hits
+  std::vector<TrackContainer> fTracks;
 
   //output file with histograms
   TString fHistoFileName {"CbmBbaAlignmentHisto.root"};
@@ -90,6 +112,11 @@ private:
   Int_t fMaxNtracks {10000};
 
   double fCostIdeal {1.e10};
+  double fCostInitial {0.};
+
+  double fChi2Total {0.};
+  long fNdfTotal {0};
+  long fFixedNdf {-1};
 
   //histograms
   TH1F* hTestHisto {nullptr};