diff --git a/macro/alignment/run_BbaAlignment_mcbm.C b/macro/alignment/run_BbaAlignment_mcbm.C
new file mode 100644
index 0000000000000000000000000000000000000000..60f8a7cb6b489f6d178b666aa394824c5cb970d3
--- /dev/null
+++ b/macro/alignment/run_BbaAlignment_mcbm.C
@@ -0,0 +1,281 @@
+/* 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();
+ // ------------------------------------------------------------------------
+}
diff --git a/reco/alignment/CMakeLists.txt b/reco/alignment/CMakeLists.txt
index 136aa977142cff387f1bd9f1a92915a6282f32c7..dcd35b5176fefa4f7f65cf5724b08e2b0214cbab 100644
--- a/reco/alignment/CMakeLists.txt
+++ b/reco/alignment/CMakeLists.txt
@@ -17,6 +17,12 @@ set(PUBLIC_DEPENDENCIES
)
set(PRIVATE_DEPENDENCIES
+ CbmMvdBase
+ CbmStsBase
+ CbmMuchBase
+ CbmTrdBase
+ CbmTofBase
+ L1
KF
bba::library
CbmStsBase
diff --git a/reco/alignment/CbmBbaAlignmentTask.cxx b/reco/alignment/CbmBbaAlignmentTask.cxx
index c8d17c92621abe2089f6c310db0bb6233984f9f9..02cb14022a225b3063ad55a3448e1f7e815b7ace 100644
--- a/reco/alignment/CbmBbaAlignmentTask.cxx
+++ b/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];
}
diff --git a/reco/alignment/CbmBbaAlignmentTask.h b/reco/alignment/CbmBbaAlignmentTask.h
index fee87f52dac7da2e4569b341c08803bea9b0217d..02bcc1cab78093c150527a0e01e4db538a1b49f3 100644
--- a/reco/alignment/CbmBbaAlignmentTask.h
+++ b/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};