BBA: a library and a dummy task for the alignment

external/CMakeLists.txt

@@ -57,6 +57,7 @@ if(DOWNLOAD_EXTERNALS)
   Add_Subdirectory(xpu)
   Set(XPU_INCLUDE_DIRECTORY ${CMAKE_CURRENT_SOURCE_DIR}/xpu/src PARENT_SCOPE)
 
+
   Include(InstallKFParticle.cmake)
   Include(InstallNicaFemto.cmake)
   Include(InstallAnalysisTree.cmake)
@@ -67,6 +68,9 @@ if(DOWNLOAD_EXTERNALS)
   Include(InstallGeometry.cmake)
 
   Include(InstallYamlCpp.cmake)
+
+  Include(InstallBBA.cmake)
+  
 else()
   # Define targets which are needed by CbmRoot but are not available
   # whithout the external packages
@@ -90,6 +94,10 @@ else()
   set_target_properties(external::fles_monitoring PROPERTIES
                         IMPORTED_LOCATION external::fles_monitoring-NOTFOUND
              )
+  add_library(bba::library STATIC IMPORTED GLOBAL)
+  set_target_properties(bba::library PROPERTIES
+                        IMPORTED_LOCATION bba::library-NOTFOUND
+            )           
 
  set(LIBRARY_OUTPUT_PATH "${CMAKE_BINARY_DIR}/lib")
  set(EXECUTABLE_OUTPUT_PATH "${CMAKE_BINARY_DIR}/bin")

external/InstallBBA.cmake

+  download_project_if_needed(PROJECT           bba
+                             GIT_REPOSITORY    "https://git.cbm.gsi.de/se.gorbunov/bba.git"
+                             GIT_TAG           "183d0e39ceeeb0ab737960a9c984192f1ba93674"
+                             SOURCE_DIR        ${CMAKE_CURRENT_SOURCE_DIR}/bba
+                             CONFIGURE_COMMAND ""
+                             BUILD_COMMAND     ""
+                             INSTALL_COMMAND   ""
+  )
+  Add_Subdirectory(bba)
+  Set(BBA_INCLUDE_DIRECTORY ${CMAKE_CURRENT_SOURCE_DIR}/bba/include PARENT_SCOPE)

macro/alignment/run_BbaAlignment.C

+/* Copyright (C) 2023-2023 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
+   SPDX-License-Identifier: GPL-3.0-only
+   Authors: S.Gorbunov[committer], N.Bluhme */
+
+// --------------------------------------------------------------------------
+//
+// Macro to run the alignment task CbmBbaAlignmentTask
+//
+// The following naming conventions are assumed:
+// Raw data file:  [dataset].raw.root
+// Transport file: [dataset].tra.root
+// Parameter file:  [dataset].par.root
+// Reconstruction file: [dataset].rec.root
+//
+// S. Gorbunov 20/01/2023
+//
+// --------------------------------------------------------------------------
+
+// commands:
+//
+// root -l -q $VMCWORKDIR/macro/alignment/run_BbaAlignment.C'("data/el.10gev.centr.eb","", "sis100_electron", "", 1)'
+//
+
+// Includes needed for IDE
+#if !defined(__CLING__)
+#include "CbmDefs.h"
+#include "CbmMCDataManager.h"
+#include "CbmMuchDigitizerQa.h"
+#include "CbmMuchHitFinderQa.h"
+#include "CbmMuchTransportQa.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(TString input = "", TString output = "", TString setup = "sis100_electron",
+                      TString paramFile = "", Int_t nEvents = -1)
+{
+
+  gROOT->SetBatch(kTRUE);
+
+  // ========================================================================
+  //          Adjust this part according to your requirements
+
+  // -----   Logger settings   ----------------------------------------------
+  FairLogger::GetLogger()->SetLogScreenLevel("INFO");
+  fair::Logger::DefineVerbosity(
+    "user1", fair::VerbositySpec::Make(fair::VerbositySpec::Info::severity, fair::VerbositySpec::Info::file_line));
+  FairLogger::GetLogger()->SetLogVerbosityLevel("user1");
+  // ------------------------------------------------------------------------
+
+  // -----   Environment   --------------------------------------------------
+  TString myName = "run_qa";                       // this macro's name for screen output
+  TString srcDir = gSystem->Getenv("VMCWORKDIR");  // top source directory
+  // ------------------------------------------------------------------------
+
+  // -----   In- and output file names   ------------------------------------
+
+  if (input.IsNull()) input = "test";
+  if (output.IsNull()) output = input;
+
+  TString traFile = input + ".tra.root";
+  TString rawFile = input + ".raw.root";
+  TString recFile = input + ".reco.root";
+
+  TString outFile = output + ".align.root";
+  TString monFile = output + ".moni_align.root";
+  if (paramFile.IsNull()) paramFile = input;
+  TString parFile = paramFile + ".par.root";
+
+  std::cout << "Inputfiles " << rawFile << " " << traFile << " " << recFile << std::endl;
+  std::cout << "Outfile " << outFile << std::endl;
+  std::cout << "Parfile " << parFile << std::endl;
+
+  // ------------------------------------------------------------------------
+
+  // -----   Load the geometry setup   -------------------------------------
+  std::cout << std::endl;
+  std::cout << "-I- " << myName << ": Loading setup " << setup << std::endl;
+  CbmSetup::Instance()->LoadSetup(setup);
+  // You can modify the pre-defined setup by using
+  // CbmSetup::Instance()->RemoveModule(ESystemId) or
+  // CbmSetup::Instance()->SetModule(ESystemId, const char*, Bool_t) or
+  // CbmSetup::Instance()->SetActive(ESystemId, Bool_t)
+  // See the class documentation of CbmSetup.
+  // ------------------------------------------------------------------------
+
+  // -----   Parameter files as input to the runtime database   -------------
+  std::cout << std::endl;
+  std::cout << "-I- " << myName << ": Defining parameter files " << 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(outFile);
+  run->SetSink(sink);
+
+  FairMonitor::GetMonitor()->EnableMonitor(kTRUE, monFile);
+  // ------------------------------------------------------------------------
+
+  // -----   MCDataManager  -----------------------------------
+  CbmMCDataManager* mcManager = new CbmMCDataManager("MCDataManager", 0);
+  mcManager->AddFile(traFile);
+
+  run->AddTask(mcManager);
+  // ------------------------------------------------------------------------
+
+  // ----- Match reco to MC ------
+  CbmMatchRecoToMC* matchTask = new CbmMatchRecoToMC();
+  run->AddTask(matchTask);
+
+  run->AddTask(new CbmTrackingDetectorInterfaceInit());  // Geometry interface initializer for tracker
+
+  // Kalman filter
+  auto kalman = new CbmKF();
+  run->AddTask(kalman);
+
+  // L1 tracking setup
+  auto l1 = new CbmL1("L1", 2, 3, 0);
+
+  // --- Material budget file names
+  TString mvdGeoTag;
+  if (CbmSetup::Instance()->GetGeoTag(ECbmModuleId::kMvd, mvdGeoTag)) {
+    TString parFile = gSystem->Getenv("VMCWORKDIR");
+    parFile += "/parameters/mvd/mvd_matbudget_" + mvdGeoTag + ".root";
+    std::cout << "Using material budget file " << parFile << std::endl;
+    l1->SetMvdMaterialBudgetFileName(parFile.Data());
+  }
+
+  TString stsGeoTag;
+  if (CbmSetup::Instance()->GetGeoTag(ECbmModuleId::kSts, stsGeoTag)) {
+    TString parFile = gSystem->Getenv("VMCWORKDIR");
+    parFile += "/parameters/sts/sts_matbudget_" + stsGeoTag + ".root";
+    std::cout << "Using material budget file " << parFile << std::endl;
+    l1->SetStsMaterialBudgetFileName(parFile.Data());
+  }
+
+  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);
+  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 " << outFile << 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/CMakeLists.txt

@@ -18,4 +18,6 @@ add_subdirectory(qa)
 add_subdirectory (tasks)
 add_subdirectory (app)
 add_subdirectory (mq)
+add_subdirectory (alignment)
+
 

reco/alignment/CMakeLists.txt

+# Create a library called "libAlignment" which includes the source files given in
+# the array .
+# The extension is already found.  Any number of sources could be listed here.
+
+set(INCLUDE_DIRECTORIES
+  ${CMAKE_CURRENT_SOURCE_DIR}
+)
+
+set(SRCS
+  CbmBbaAlignmentTask.cxx
+)
+
+set(LIBRARY_NAME CbmAlignment)
+set(LINKDEF ${LIBRARY_NAME}LinkDef.h)
+set(PUBLIC_DEPENDENCIES
+  FairRoot::Base
+)
+
+set(PRIVATE_DEPENDENCIES
+  L1
+  bba::library
+)
+
+generate_cbm_library()

reco/alignment/CbmAlignmentLinkDef.h

+/* Copyright (C) 2006-2020 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
+   SPDX-License-Identifier: GPL-3.0-only
+   Authors: Denis Bertini [committer], Volker Friese */
+
+#ifdef __CINT__
+
+#pragma link off all globals;
+#pragma link off all classes;
+#pragma link off all functions;
+
+
+#pragma link C++ class CbmBbaAlignmentTask + ;
+
+
+#endif

reco/alignment/CbmBbaAlignmentTask.cxx

+/* Copyright (C) 2023-2023 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
+   SPDX-License-Identifier: GPL-3.0-only
+   Authors: S.Gorbunov[committer], N.Bluhme */
+
+/// @file    CbmBbaAlignmentTask.cxx
+/// @author  Sergey Gorbunov
+/// @author  Nora Bluhme
+/// @date    20.01.2023
+/// @brief   Task class for alignment
+///
+
+// Cbm Headers ----------------------
+#include "CbmBbaAlignmentTask.h"
+
+#include "CbmL1PFFitter.h"
+
+// ROOT headers
+
+#include "FairRootManager.h"
+
+#include "TClonesArray.h"
+#include "TFile.h"
+#include "TH1F.h"
+
+//
+
+#include "CbmStsTrack.h"
+
+#include "bba/BBA.h"
+
+// c++ and std headers
+
+#include <iostream>
+
+
+CbmBbaAlignmentTask::CbmBbaAlignmentTask(const char* name, Int_t iVerbose, TString histoFileName)
+  : FairTask(name, iVerbose)
+  , fHistoFileName(histoFileName)
+{
+  TFile* curFile           = gFile;
+  TDirectory* curDirectory = gDirectory;
+
+  if (!(fHistoFileName == "")) { fHistoFile = new TFile(fHistoFileName.Data(), "RECREATE"); }
+  else {
+    fHistoFile = gFile;
+  }
+
+  fHistoFile->cd();
+
+  fHistoDir = fHistoFile->mkdir("Bba");
+  fHistoDir->cd();
+
+  hTestHisto = new TH1F("hTestHisto", "hTestHisto", 100, 0., 1.);
+
+  gFile      = curFile;
+  gDirectory = curDirectory;
+}
+
+
+CbmBbaAlignmentTask::~CbmBbaAlignmentTask() {}
+
+
+InitStatus CbmBbaAlignmentTask::Init()
+{
+  //Get ROOT Manager
+  FairRootManager* ioman = FairRootManager::Instance();
+
+  if (!ioman) {
+    LOG(error) << "CbmBbaAlignmentTask::Init :: RootManager not instantiated!";
+    return kERROR;
+  }
+
+  // Get sts tracks
+  fStsTrackArray = (TClonesArray*) ioman->GetObject("StsTrack");
+  if (!fStsTrackArray) {
+    LOG(error) << "CbmBbaAlignmentTask::Init: Sts track-array not found!";
+    return kERROR;
+  }
+
+  // MC track match
+  fMCTrackArray = (TClonesArray*) ioman->GetObject("MCTrack");
+  if (!fMCTrackArray) {
+    Warning("CbmBbaAlignmentTask::Init", "mc track array not found!");
+    return kERROR;
+  }
+
+  // Track match
+  fStsTrackMatchArray = (TClonesArray*) ioman->GetObject("StsTrackMatch");
+  if (fStsTrackMatchArray == 0) {
+    LOG(error) << "CbmBbaAlignmentTask::Init: track match array not found!";
+    return kERROR;
+  }
+
+  return kSUCCESS;
+}
+
+void CbmBbaAlignmentTask::Exec(Option_t* /*opt*/)
+{
+
+  std::cout << "BBA: exec event N " << fNEvents << std::endl;
+
+  fNEvents++;
+
+  bba::BBA alignment;
+
+  alignment.printInfo();
+
+  // Check fit quality of the STS tracks
+
+  std::vector<CbmStsTrack> vTracks;
+  vTracks.reserve(fStsTrackArray->GetEntriesFast());
+
+  for (int iTr = 0; iTr < fStsTrackArray->GetEntriesFast(); iTr++) {
+    CbmStsTrack* stsTrack = ((CbmStsTrack*) fStsTrackArray->At(iTr));
+    const auto* par       = stsTrack->GetParamFirst();
+    if (!std::isfinite(par->GetQp())) continue;
+    if (fabs(par->GetQp()) > 1.) continue;
+    vTracks.push_back(*stsTrack);
+  }
+
+  std::vector<int> pdg(vTracks.size(), 211);  // pion PDG
+
+  CbmL1PFFitter fitter;
+  fitter.Fit(vTracks, pdg);
+
+  double chi2Total  = 0;
+  long int ndfTotal = 1;
+
+  for (unsigned int iTr = 0; iTr < vTracks.size(); iTr++) {
+    if (!std::isfinite(vTracks[iTr].GetChiSq())) continue;
+    chi2Total += vTracks[iTr].GetChiSq();
+    ndfTotal += vTracks[iTr].GetNDF();
+  }
+  std::cout << "BBA: chi2/ndf = " << chi2Total / ndfTotal << std::endl;
+}
+
+void CbmBbaAlignmentTask::Finish()
+{
+  TDirectory* curr   = gDirectory;
+  TFile* currentFile = gFile;
+  // Open output file and write histograms
+
+  fHistoFile->cd();
+  WriteHistosCurFile(fHistoDir);
+  if (!(fHistoFileName == "")) {
+    fHistoFile->Close();
+    fHistoFile->Delete();
+  }
+  gFile      = currentFile;
+  gDirectory = curr;
+}
+
+void CbmBbaAlignmentTask::WriteHistosCurFile(TObject* obj)
+{
+  if (!obj->IsFolder()) obj->Write();
+  else {
+    TDirectory* cur    = gDirectory;
+    TFile* currentFile = gFile;
+
+    TDirectory* sub = cur->GetDirectory(obj->GetName());
+    sub->cd();
+    TList* listSub = (static_cast<TDirectory*>(obj))->GetList();
+    TIter it(listSub);
+    while (TObject* obj1 = it())
+      WriteHistosCurFile(obj1);
+    cur->cd();
+    gFile      = currentFile;
+    gDirectory = cur;
+  }
+}
+
+ClassImp(CbmBbaAlignmentTask);

reco/alignment/CbmBbaAlignmentTask.h

+/* Copyright (C) 2023-2023 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
+   SPDX-License-Identifier: GPL-3.0-only
+   Authors: S.Gorbunov[committer], N.Bluhme */
+
+/// @file    CbmBbaAlignmentTask.h
+/// @author  Sergey Gorbunov
+/// @author  Nora Bluhme
+/// @date    20.01.2023
+/// @brief   Task class for alignment
+///
+
+#ifndef CbmBbaAlignmentTask_H
+#define CbmBbaAlignmentTask_H
+
+#include "FairTask.h"
+
+class TClonesArray;
+class TFile;
+class TDirectory;
+class TH1F;
+
+class CbmBbaAlignmentTask : public FairTask {
+public:
+  // Constructors/Destructors ---------
+  CbmBbaAlignmentTask(const char* name = "CbmBbaAlignmentTask", Int_t iVerbose = 0,
+                      TString histoFileName = "CbmBbaAlignmentHisto.root");
+  ~CbmBbaAlignmentTask();
+
+  Int_t GetZtoNStation(Double_t getZ);
+
+  InitStatus Init();
+  void Exec(Option_t* opt);
+  void Finish();
+
+private:
+  const CbmBbaAlignmentTask& operator=(const CbmBbaAlignmentTask&);
+  CbmBbaAlignmentTask(const CbmBbaAlignmentTask&);
+
+  void WriteHistosCurFile(TObject* obj);
+  int GetHistoIndex(int pdg);
+
+  //input branches
+  TClonesArray* fStsTrackArray {nullptr};
+  TClonesArray* fMCTrackArray {nullptr};
+  TClonesArray* fStsTrackMatchArray {nullptr};
+
+  //output file with histograms
+  TString fHistoFileName {"CbmBbaAlignmentHisto.root"};
+  TFile* fHistoFile {nullptr};
+  TDirectory* fHistoDir {nullptr};
+
+  Int_t fNEvents {0};
+
+  //histograms
+  TH1F* hTestHisto {nullptr};
+
+  ClassDef(CbmBbaAlignmentTask, 1);
+};
+
+#endif