From e9f6bf899a2ee434b69a77f292aa7eb2bcaa8797 Mon Sep 17 00:00:00 2001
From: "se.gorbunov" <se.gorbunov@gsi.de>
Date: Tue, 28 Feb 2023 19:17:46 +0000
Subject: [PATCH] bba: draft version of the alignmnet task
---
reco/L1/ParticleFinder/CbmL1PFFitter.cxx | 56 +++--
reco/L1/ParticleFinder/CbmL1PFFitter.h | 6 +-
reco/alignment/CMakeLists.txt | 1 +
reco/alignment/CbmBbaAlignmentTask.cxx | 287 +++++++++++++++++++++--
reco/alignment/CbmBbaAlignmentTask.h | 42 +++-
5 files changed, 347 insertions(+), 45 deletions(-)
diff --git a/reco/L1/ParticleFinder/CbmL1PFFitter.cxx b/reco/L1/ParticleFinder/CbmL1PFFitter.cxx
index 3264abf8aa..f95f3c6698 100644
--- a/reco/L1/ParticleFinder/CbmL1PFFitter.cxx
+++ b/reco/L1/ParticleFinder/CbmL1PFFitter.cxx
@@ -46,6 +46,7 @@
//typedef L1Fit1 L1Fit;
using std::vector;
+using namespace std;
namespace NS_L1TrackFitter
{
@@ -100,7 +101,7 @@ inline CbmL1PFFitter::PFFieldRegion::PFFieldRegion(const L1FieldRegion& fld, int
void FilterFirst(L1Fit& fit, fvec& x, fvec& y, fvec& dxx, fvec& dxy, fvec& dyy)
{
L1TrackPar& tr = fit.Tr();
- tr.ResetErrors(dxx, dyy, vINF, vINF, 1., 1.e6, 1.e2);
+ tr.ResetErrors(dxx, dyy, 1., 1., 1., 1.e6, 1.e2);
tr.C10 = dxy;
tr.x = x;
tr.y = y;
@@ -109,18 +110,13 @@ void FilterFirst(L1Fit& fit, fvec& x, fvec& y, fvec& dxx, fvec& dxy, fvec& dyy)
tr.NDF = -3.0;
}
-void CbmL1PFFitter::Fit(vector<CbmStsTrack>& Tracks, vector<int>& pidHypo)
+void CbmL1PFFitter::Fit(std::vector<CbmStsTrack>& Tracks, const std::vector<CbmMvdHit>& vMvdHits,
+ const std::vector<CbmStsHit>& vStsHits, const std::vector<int>& pidHypo)
{
L1FieldValue b0, b1, b2 _fvecalignment;
L1FieldRegion fld _fvecalignment;
-
-
- FairRootManager* fManger = FairRootManager::Instance();
- TClonesArray* listStsHits = (TClonesArray*) fManger->GetObject("StsHit");
- int NMvdStations = CbmL1::Instance()->fpAlgo->GetNstationsBeforePipe();
- TClonesArray* listMvdHits = 0;
- if (NMvdStations > 0.) listMvdHits = (TClonesArray*) fManger->GetObject("MvdHit");
+ // fld.SetUseOriginalField();
static int nHits = CbmL1::Instance()->fpAlgo->GetParameters()->GetNstationsActive();
int iVec = 0, i = 0;
@@ -169,7 +165,7 @@ void CbmL1PFFitter::Fit(vector<CbmStsTrack>& Tracks, vector<int>& pidHypo)
T.t = 0.;
T.vi = 0.;
- T.ResetErrors(1.e2, 1.e2, 1.e4, 1.e4, 1.e4, 1.e6, 1.e2);
+ T.ResetErrors(1.e2, 1.e2, 1., 1., 1., 1.e6, 1.e2);
if (N_vTracks - itrack < static_cast<unsigned short>(fvec::size())) nTracks_SIMD = N_vTracks - itrack;
for (i = 0; i < nTracks_SIMD; i++) {
@@ -179,6 +175,8 @@ void CbmL1PFFitter::Fit(vector<CbmStsTrack>& Tracks, vector<int>& pidHypo)
T.tx[i] = tr[i]->GetParamFirst()->GetTx();
T.ty[i] = tr[i]->GetParamFirst()->GetTy();
T.qp[i] = tr[i]->GetParamFirst()->GetQp();
+ T.t[i] = 0.;
+ T.vi[i] = 0.;
T.z[i] = tr[i]->GetParamFirst()->GetZ();
T.C00[i] = tr[i]->GetParamFirst()->GetCovariance(0, 0);
T.C10[i] = tr[i]->GetParamFirst()->GetCovariance(1, 0);
@@ -218,9 +216,8 @@ void CbmL1PFFitter::Fit(vector<CbmStsTrack>& Tracks, vector<int>& pidHypo)
float posx = 0.f, posy = 0.f, posz = 0.f, time = 0.;
if (i < nHitsTrackMvd) {
- if (!listMvdHits) continue;
- int hitIndex = tr[iVec]->GetMvdHitIndex(i);
- CbmMvdHit* hit = L1_DYNAMIC_CAST<CbmMvdHit*>(listMvdHits->At(hitIndex));
+ int hitIndex = tr[iVec]->GetMvdHitIndex(i);
+ const CbmMvdHit* hit = &(vMvdHits[hitIndex]);
posx = hit->GetX();
posy = hit->GetY();
@@ -236,9 +233,8 @@ void CbmL1PFFitter::Fit(vector<CbmStsTrack>& Tracks, vector<int>& pidHypo)
dt2[ista][iVec] = hit->GetTimeError() * hit->GetTimeError();
}
else {
- if (!listStsHits) continue;
- int hitIndex = tr[iVec]->GetHitIndex(i - nHitsTrackMvd);
- CbmStsHit* hit = L1_DYNAMIC_CAST<CbmStsHit*>(listStsHits->At(hitIndex));
+ int hitIndex = tr[iVec]->GetHitIndex(i - nHitsTrackMvd);
+ const CbmStsHit* hit = &(vStsHits[hitIndex]);
posx = hit->GetX();
posy = hit->GetY();
@@ -291,7 +287,9 @@ void CbmL1PFFitter::Fit(vector<CbmStsTrack>& Tracks, vector<int>& pidHypo)
}
}
}
+
// fit forward
+
i = 0;
FilterFirst(fit, x_first, y_first, fstC00, fstC10, fstC11);
@@ -359,7 +357,8 @@ void CbmL1PFFitter::Fit(vector<CbmStsTrack>& Tracks, vector<int>& pidHypo)
tr[iVec]->SetParamLast(&par);
}
- //fit backward
+ // fit backward
+
fit.SetQp0(T.qp);
i = nHits - 1;
@@ -441,6 +440,29 @@ void CbmL1PFFitter::Fit(vector<CbmStsTrack>& Tracks, vector<int>& pidHypo)
delete[] fB;
}
+void CbmL1PFFitter::Fit(vector<CbmStsTrack>& Tracks, const vector<int>& pidHypo)
+{
+
+ FairRootManager* fManger = FairRootManager::Instance();
+ TClonesArray* mvdHitArray = (TClonesArray*) fManger->GetObject("MvdHit");
+ TClonesArray* stsHitArray = (TClonesArray*) fManger->GetObject("StsHit");
+
+ std::vector<CbmMvdHit> vMvdHits;
+ std::vector<CbmStsHit> vStsHits;
+
+ for (int ih = 0; ih < mvdHitArray->GetEntriesFast(); ih++) {
+ CbmMvdHit hit = *dynamic_cast<const CbmMvdHit*>(mvdHitArray->At(ih));
+ vMvdHits.push_back(hit);
+ }
+
+ for (int ih = 0; ih < stsHitArray->GetEntriesFast(); ih++) {
+ vStsHits.push_back(*dynamic_cast<const CbmStsHit*>(stsHitArray->At(ih)));
+ }
+
+ Fit(Tracks, vMvdHits, vStsHits, pidHypo);
+}
+
+
void CbmL1PFFitter::GetChiToVertex(vector<CbmStsTrack>& Tracks, vector<PFFieldRegion>& field, vector<float>& chiToVtx,
CbmKFVertex& primVtx, float chiPrim)
{
diff --git a/reco/L1/ParticleFinder/CbmL1PFFitter.h b/reco/L1/ParticleFinder/CbmL1PFFitter.h
index 2189862fad..1955f4db34 100644
--- a/reco/L1/ParticleFinder/CbmL1PFFitter.h
+++ b/reco/L1/ParticleFinder/CbmL1PFFitter.h
@@ -22,6 +22,8 @@
#include <vector>
+class CbmMvdHit;
+class CbmStsHit;
class CbmStsTrack;
class L1FieldRegion;
class CbmKFVertex;
@@ -42,7 +44,9 @@ public:
~CbmL1PFFitter();
//functions for fitting CbmStsTrack
- void Fit(std::vector<CbmStsTrack>& Tracks, std::vector<int>& pidHypo);
+ void Fit(std::vector<CbmStsTrack>& Tracks, const std::vector<CbmMvdHit>& vMvdHits,
+ const std::vector<CbmStsHit>& vStsHits, const std::vector<int>& pidHypo);
+ void Fit(std::vector<CbmStsTrack>& Tracks, const std::vector<int>& pidHypo);
void CalculateFieldRegion(std::vector<CbmStsTrack>& Tracks, std::vector<PFFieldRegion>& Field);
void CalculateFieldRegionAtLastPoint(std::vector<CbmStsTrack>& Tracks, std::vector<PFFieldRegion>& field);
void GetChiToVertex(std::vector<CbmStsTrack>& Tracks, std::vector<PFFieldRegion>& field, std::vector<float>& chiToVtx,
diff --git a/reco/alignment/CMakeLists.txt b/reco/alignment/CMakeLists.txt
index 196318ae96..c2e9448ef0 100644
--- a/reco/alignment/CMakeLists.txt
+++ b/reco/alignment/CMakeLists.txt
@@ -19,6 +19,7 @@ set(PUBLIC_DEPENDENCIES
set(PRIVATE_DEPENDENCIES
L1
bba::library
+ CbmStsBase
)
generate_cbm_library()
diff --git a/reco/alignment/CbmBbaAlignmentTask.cxx b/reco/alignment/CbmBbaAlignmentTask.cxx
index c2cee83ad8..9354603a7f 100644
--- a/reco/alignment/CbmBbaAlignmentTask.cxx
+++ b/reco/alignment/CbmBbaAlignmentTask.cxx
@@ -21,9 +21,13 @@
#include "TClonesArray.h"
#include "TFile.h"
#include "TH1F.h"
+#include "TRandom.h"
//
+#include "CbmMvdHit.h"
+#include "CbmStsHit.h"
+#include "CbmStsSetup.h"
#include "CbmStsTrack.h"
#include "bba/BBA.h"
@@ -70,27 +74,36 @@ InitStatus CbmBbaAlignmentTask::Init()
return kERROR;
}
+ // Get hits
+
+ fInputMvdHits = (TClonesArray*) ioman->GetObject("MvdHit");
+ fInputStsHits = (TClonesArray*) ioman->GetObject("StsHit");
+
// Get sts tracks
- fStsTrackArray = (TClonesArray*) ioman->GetObject("StsTrack");
- if (!fStsTrackArray) {
+ fInputStsTracks = (TClonesArray*) ioman->GetObject("StsTrack");
+ if (!fInputStsTracks) {
LOG(error) << "CbmBbaAlignmentTask::Init: Sts track-array not found!";
return kERROR;
}
// MC track match
- fMCTrackArray = (TClonesArray*) ioman->GetObject("MCTrack");
- if (!fMCTrackArray) {
+ fInputMcTracks = (TClonesArray*) ioman->GetObject("MCTrack");
+ if (!fInputMcTracks) {
Warning("CbmBbaAlignmentTask::Init", "mc track array not found!");
return kERROR;
}
// Track match
- fStsTrackMatchArray = (TClonesArray*) ioman->GetObject("StsTrackMatch");
- if (fStsTrackMatchArray == 0) {
+ fInputStsTrackMatches = (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;
}
@@ -101,41 +114,269 @@ void CbmBbaAlignmentTask::Exec(Option_t* /*opt*/)
fNEvents++;
- bba::BBA alignment;
+ if ((int) fTracks.size() >= fMaxNtracks) { return; }
- alignment.printInfo();
+ // select STS tracks for alignment and store them
- // Check fit quality of the STS tracks
+ for (int iTr = 0; iTr < fInputStsTracks->GetEntriesFast(); iTr++) {
- std::vector<CbmStsTrack> vTracks;
- vTracks.reserve(fStsTrackArray->GetEntriesFast());
+ if ((int) fTracks.size() >= fMaxNtracks) { break; }
- for (int iTr = 0; iTr < fStsTrackArray->GetEntriesFast(); iTr++) {
- CbmStsTrack* stsTrack = ((CbmStsTrack*) fStsTrackArray->At(iTr));
- const auto* par = stsTrack->GetParamFirst();
+ CbmStsTrack* stsTrack = ((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;
- vTracks.push_back(*stsTrack);
+ 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);
+ }
+
+ 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);
+ }
+
+ fTracks.push_back(track);
+
+ } // tracks
+}
+
+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());
+
+ assert(iStation >= 0 && iStation < CbmStsSetup::Instance()->GetNofStations());
+
+ double dx = par[3 * iStation + 0];
+ double dy = par[3 * iStation + 1];
+ double dz = par[3 * iStation + 2];
+
+ CbmStsHit& hitAligned = fStsHitsAligned[ih];
+
+ hitAligned.SetX(hitOrig.GetX() + dx);
+ hitAligned.SetY(hitOrig.GetY() + dy);
+ hitAligned.SetZ(hitOrig.GetZ() + dz);
}
+}
+
+double CbmBbaAlignmentTask::CostFunction(const std::vector<double>& par)
+{
+
+ // apply new parameters to the hits
+
+ ApplyAlignment(par);
- std::vector<int> pdg(vTracks.size(), 211); // pion PDG
+ auto newTracks = fTracks;
CbmL1PFFitter fitter;
- fitter.Fit(vTracks, pdg);
+ fitter.Fit(newTracks, fMvdHits, fStsHitsAligned, fTrackPdg);
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();
+ 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();
}
- std::cout << "BBA: chi2/ndf = " << chi2Total / ndfTotal << std::endl;
+
+ double cost = chi2Total / ndfTotal;
+ std::cout << "BBA: cost function: n tracks " << nGoodTracks << ", cost " << cost
+ << ", diff to ideal cost: " << cost - fCostIdeal << std::endl;
+ return cost;
+ if (nGoodTracks == (int) fTracks.size()) { return cost; }
+ return 1.e30;
}
void CbmBbaAlignmentTask::Finish()
{
+ //
+ // perform the alignment
+ //
+
+ std::cout << "BBA: start the alignment procedure with " << fTracks.size() << " tracks ..." << std::endl;
+
+ // init auxiliary arrays
+
+ fStsHitsAligned = fStsHits;
+ fTrackPdg.resize(fTracks.size(), 211); // pion PDG
+
+ // init alignmemt module
+
+ int nStsStations = CbmStsSetup::Instance()->GetNofStations();
+ int nParameters = 3 * nStsStations; // x,yz
+
+ std::vector<bba::Parameter> par(nParameters);
+
+ for (int ip = 0; ip < nParameters; ip++) {
+ bba::Parameter& p = par[ip];
+ p.SetActive(0);
+ p.SetValue(0.);
+ p.SetBoundaryMin(-2.); // +-1 cm alignment range
+ p.SetBoundaryMax(2.);
+ p.SetStepMin(1.e-6);
+ p.SetStepMax(0.1);
+ 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[3 * 4 + 0].SetActive(1);
+ par[3 * 4 + 1].SetActive(1);
+ par[3 * 4 + 2].SetActive(1);
+
+ 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);
+
+
+ gRandom->SetSeed(1);
+
+ for (int is = 0; is < nStsStations; is++) {
+ bba::Parameter& px = par[3 * is + 0];
+ bba::Parameter& py = par[3 * is + 1];
+ bba::Parameter& pz = par[3 * is + 2];
+ //py.SetActive(1);
+ //pz.SetStepMin(1.e-4);
+ // +- 0.5 cm misalignment
+
+ //if (px.IsActive()) { px.SetValue(.1); } //gRandom->Uniform(2 * d) - d); }
+ //if (py.IsActive()) { py.SetValue(.1); } //gRandom->Uniform(2. * d) - d); }
+ //if (pz.IsActive()) { pz.SetValue(.1); } //gRandom->Uniform(2. * d) - d); }
+ double d = 0.5;
+ if (px.IsActive()) { px.SetValue(gRandom->Uniform(2. * d) - d); }
+ if (py.IsActive()) { py.SetValue(gRandom->Uniform(2. * d) - d); }
+ if (pz.IsActive()) { pz.SetValue(gRandom->Uniform(2. * d) - d); }
+ }
+
+ bba::BBA alignment;
+
+ alignment.setParameters(par);
+
+ auto lambdaCostFunction = [this](const std::vector<double>& p) { return CostFunction(p); };
+
+ alignment.setChi2(lambdaCostFunction);
+ alignment.printInfo();
+
+ std::vector<double> parAligned;
+ std::vector<double> parMisaligned;
+ {
+ const std::vector<double>& r = alignment.getResult();
+ for (unsigned int i = 0; i < r.size(); i++) {
+ parMisaligned.push_back(r[i]);
+ parAligned.push_back(0.);
+ }
+ }
+
+ {
+ 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;
+ 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();
+ }
+
+ std::cout << "Initial nTracks " << nGoodTracks << " chi2/ndf " << chi2Total / ndfTotal << std::endl;
+ }
+
+ fCostIdeal = CostFunction(parAligned);
+
+ std::cout << " cost function for the true parameters: " << fCostIdeal << std::endl;
+
+ alignment.align();
+
+ std::cout << " cost function for the true parameters: " << fCostIdeal << std::endl;
+
+ std::cout << " Misaligned parameters: " << std::endl;
+ for (int is = 0; is < nStsStations; is++) {
+ const std::vector<double>& r = parMisaligned;
+ std::cout << "Sts station " << is << ": x " << r[3 * is + 0] << " y " << r[3 * is + 1] << " z " << r[3 * is + 2]
+ << std::endl;
+ }
+
+ std::cout << " Alignment results: " << std::endl;
+
+ for (int is = 0; is < nStsStations; is++) {
+ const std::vector<double>& r = alignment.getResult();
+ std::cout << "Sts station " << is << ": x " << r[3 * is + 0] << " y " << r[3 * is + 1] << " z " << r[3 * is + 2]
+ << std::endl;
+ }
+
+ // store the histograms
+
TDirectory* curr = gDirectory;
TFile* currentFile = gFile;
// Open output file and write histograms
diff --git a/reco/alignment/CbmBbaAlignmentTask.h b/reco/alignment/CbmBbaAlignmentTask.h
index a316f334fb..fb24957a6b 100644
--- a/reco/alignment/CbmBbaAlignmentTask.h
+++ b/reco/alignment/CbmBbaAlignmentTask.h
@@ -12,6 +12,10 @@
#ifndef CbmBbaAlignmentTask_H
#define CbmBbaAlignmentTask_H
+#include "CbmMvdHit.h"
+#include "CbmStsHit.h"
+#include "CbmStsTrack.h"
+
#include "FairTask.h"
class TClonesArray;
@@ -19,6 +23,14 @@ class TFile;
class TDirectory;
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 ---------
@@ -39,10 +51,28 @@ private:
void WriteHistosCurFile(TObject* obj);
int GetHistoIndex(int pdg);
- //input branches
- TClonesArray* fStsTrackArray {nullptr};
- TClonesArray* fMCTrackArray {nullptr};
- TClonesArray* fStsTrackMatchArray {nullptr};
+ 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};
+
+ TClonesArray* fInputMcTracks {nullptr}; // Mc info for debugging
+ TClonesArray* fInputStsTrackMatches {nullptr}; // Mc info for debugging
+
+ // collection of selected tracks and hits
+ std::vector<CbmStsTrack> fTracks;
+ std::vector<CbmMvdHit> fMvdHits;
+ std::vector<CbmStsHit> fStsHits;
+
+ // temporary array with aligned hits
+ std::vector<CbmStsHit> fStsHitsAligned;
+
+ // array with pdg hypothesis for tracks
+ std::vector<int> fTrackPdg;
//output file with histograms
TString fHistoFileName {"CbmBbaAlignmentHisto.root"};
@@ -51,6 +81,10 @@ private:
Int_t fNEvents {0};
+ Int_t fMaxNtracks {10000};
+
+ double fCostIdeal {1.e10};
+
//histograms
TH1F* hTestHisto {nullptr};
--
GitLab