Skip to content
Snippets Groups Projects
Commit 2ba38019 authored by Sergey Gorbunov's avatar Sergey Gorbunov Committed by Sergey Gorbunov
Browse files

BBA: multithreading, fix full-CBM mode

parent 3fd35831
No related branches found
No related tags found
1 merge request!1682BBA: multithreading, fix full-CBM mode
Pipeline #27500 passed
Stage: build
    Stage: package
      Stage: verify
        Hide whitespace changes
        Inline Side-by-side
        macro/alignment/run_BbaAlignment.C
        + 3
        − 3
        View file @ 2ba38019
        ...@@ -149,9 +149,9 @@ void run_BbaAlignment(TString input = "", TString output = "", TString setup = " ...@@ -149,9 +149,9 @@ void run_BbaAlignment(TString input = "", TString output = "", TString setup = "
        run->AddTask(l1); run->AddTask(l1);
        // ----- BBA alignment -------------------------------------------- // ----- BBA alignment --------------------------------------------
        CbmBbaAlignmentTask* alignemnt = new CbmBbaAlignmentTask(); CbmBbaAlignmentTask* alignment = new CbmBbaAlignmentTask();
        //CbmKFTrackQa* alignemnt = new CbmKFTrackQa(); alignment->SetStsTrackingMode();
        run->AddTask(alignemnt); run->AddTask(alignment);
        // ----- Parameter database -------------------------------------------- // ----- Parameter database --------------------------------------------
        ......
        macro/alignment/run_BbaAlignment_mcbm.C
        + 5
        − 5
        View file @ 2ba38019
        ...@@ -103,7 +103,7 @@ void run_BbaAlignment_mcbm(Int_t nEvents = 0, TString dataset = "data/mcbm_beam_ ...@@ -103,7 +103,7 @@ void run_BbaAlignment_mcbm(Int_t nEvents = 0, TString dataset = "data/mcbm_beam_
        TString geoTag; TString geoTag;
        // - MUCH digitisation parameters // - MUCH digitisation parameters
        TString muchParFile {}; TString muchParFile{};
        if (setup->GetGeoTag(ECbmModuleId::kMuch, geoTag)) { if (setup->GetGeoTag(ECbmModuleId::kMuch, geoTag)) {
        bool mcbmFlag = geoTag.Contains("mcbm", TString::kIgnoreCase); bool mcbmFlag = geoTag.Contains("mcbm", TString::kIgnoreCase);
        muchParFile = srcDir + "/parameters/much/much_"; muchParFile = srcDir + "/parameters/much/much_";
        ...@@ -161,7 +161,7 @@ void run_BbaAlignment_mcbm(Int_t nEvents = 0, TString dataset = "data/mcbm_beam_ ...@@ -161,7 +161,7 @@ void run_BbaAlignment_mcbm(Int_t nEvents = 0, TString dataset = "data/mcbm_beam_
        FairRootFileSink* sink = new FairRootFileSink(sinkFile); FairRootFileSink* sink = new FairRootFileSink(sinkFile);
        run->SetSink(sink); run->SetSink(sink);
        TString monitorFile {sinkFile}; TString monitorFile{sinkFile};
        monitorFile.ReplaceAll("qa", "qa.monitor"); monitorFile.ReplaceAll("qa", "qa.monitor");
        FairMonitor::GetMonitor()->EnableMonitor(kTRUE, monitorFile); FairMonitor::GetMonitor()->EnableMonitor(kTRUE, monitorFile);
        // ------------------------------------------------------------------------ // ------------------------------------------------------------------------
        ...@@ -210,9 +210,9 @@ void run_BbaAlignment_mcbm(Int_t nEvents = 0, TString dataset = "data/mcbm_beam_ ...@@ -210,9 +210,9 @@ void run_BbaAlignment_mcbm(Int_t nEvents = 0, TString dataset = "data/mcbm_beam_
        run->AddTask(l1); run->AddTask(l1);
        // ----- BBA alignment -------------------------------------------- // ----- BBA alignment --------------------------------------------
        CbmBbaAlignmentTask* alignemnt = new CbmBbaAlignmentTask(); CbmBbaAlignmentTask* alignment = new CbmBbaAlignmentTask();
        //CbmKFTrackQa* alignemnt = new CbmKFTrackQa(); alignment->SetMcbmTrackingMode();
        run->AddTask(alignemnt); run->AddTask(alignment);
        // ----- Parameter database -------------------------------------------- // ----- Parameter database --------------------------------------------
        ......
        reco/KF/CbmKfTrackFitter.cxx
        + 1
        − 2
        View file @ 2ba38019
        ...@@ -53,7 +53,7 @@ void CbmKfTrackFitter::Track::MakeConsistent() ...@@ -53,7 +53,7 @@ void CbmKfTrackFitter::Track::MakeConsistent()
        fFirstHitNode = fNodes.size() - 1; fFirstHitNode = fNodes.size() - 1;
        fLastHitNode = 0; fLastHitNode = 0;
        for (int i = 0; i < (int) fNodes.size(); i++) { for (int i = 0; i < (int) fNodes.size(); i++) {
        if (fNodes[i].fMxy.NdfX()[0] + fNodes[i].fMxy.NdfY()[0] > 0) { if (fNodes[i].fIsXySet) {
        fFirstHitNode = std::min(fFirstHitNode, i); fFirstHitNode = std::min(fFirstHitNode, i);
        fLastHitNode = std::max(fLastHitNode, i); fLastHitNode = std::max(fLastHitNode, i);
        } }
        ...@@ -488,7 +488,6 @@ void CbmKfTrackFitter::FitTrack(CbmKfTrackFitter::Track& t) ...@@ -488,7 +488,6 @@ void CbmKfTrackFitter::FitTrack(CbmKfTrackFitter::Track& t)
        fFit.SetMask(fmask::One()); fFit.SetMask(fmask::One());
        fFit.SetParticleMass(fMass); fFit.SetParticleMass(fMass);
        t.fIsMsQp0Set = false;
        ca::FieldRegion<ca::fvec> field _fvecalignment; ca::FieldRegion<ca::fvec> field _fvecalignment;
        field.SetUseOriginalField(); field.SetUseOriginalField();
        ......
        reco/L1/catools/CaToolsField.h
        + 7
        − 0
        View file @ 2ba38019
        ...@@ -13,6 +13,8 @@ ...@@ -13,6 +13,8 @@
        #include "FairField.h" #include "FairField.h"
        #include "FairRunAna.h" #include "FairRunAna.h"
        #include <mutex>
        namespace cbm::ca::tools namespace cbm::ca::tools
        { {
        /// pass the original magnetic field to L1Algo /// pass the original magnetic field to L1Algo
        ...@@ -23,7 +25,12 @@ namespace cbm::ca::tools ...@@ -23,7 +25,12 @@ namespace cbm::ca::tools
        assert(FairRunAna::Instance()->GetField()); assert(FairRunAna::Instance()->GetField());
        double pos[3] = {x, y, z}; double pos[3] = {x, y, z};
        double B[3] = {0., 0., 0.}; double B[3] = {0., 0., 0.};
        // protect the field access
        // TODO: make CbmField thread-safe
        static std::mutex mymutex;
        mymutex.lock();
        FairRunAna::Instance()->GetField()->GetFieldValue(pos, B); FairRunAna::Instance()->GetField()->GetFieldValue(pos, B);
        mymutex.unlock();
        return std::tuple(B[0], B[1], B[2]); return std::tuple(B[0], B[1], B[2]);
        }; };
        ca::FieldRegion<ca::fvec>::SetOriginalField(fld); ca::FieldRegion<ca::fvec>::SetOriginalField(fld);
        ......
        reco/alignment/CbmBbaAlignmentTask.cxx
        + 142
        − 65
        View file @ 2ba38019
        ...@@ -41,10 +41,8 @@ ...@@ -41,10 +41,8 @@
        #include "bba/BBA.h" #include "bba/BBA.h"
        #include <cmath> #include <cmath>
        // c++ and std headers
        #include <iostream> #include <iostream>
        #include <thread>
        namespace namespace
        { {
        ...@@ -101,6 +99,22 @@ CbmBbaAlignmentTask::~CbmBbaAlignmentTask() {} ...@@ -101,6 +99,22 @@ CbmBbaAlignmentTask::~CbmBbaAlignmentTask() {}
        InitStatus CbmBbaAlignmentTask::Init() InitStatus CbmBbaAlignmentTask::Init()
        { {
        {
        const char* env = std::getenv("OMP_NUM_THREADS");
        if (env) {
        fNthreads = TString(env).Atoi();
        LOG(info) << " found environment variable OMP_NUM_THREADS = \"" << env << "\", read as integer: " << fNthreads;
        }
        }
        // ensure that at least one thread is set
        if (fNthreads < 1) {
        fNthreads = 1;
        }
        //fNthreads = 1;
        fFitter.Init(); fFitter.Init();
        //Get ROOT Manager //Get ROOT Manager
        ...@@ -170,14 +184,26 @@ InitStatus CbmBbaAlignmentTask::Init() ...@@ -170,14 +184,26 @@ InitStatus CbmBbaAlignmentTask::Init()
        fHistoDir->mkdir(n1); fHistoDir->mkdir(n1);
        fHistoDir->cd(n1); fHistoDir->cd(n1);
        double sizeX = 0.1;
        double sizeY = 0.1;
        if (fTrackingMode == kSts) {
        sizeX = 0.1;
        sizeY = 0.1;
        }
        else if (fTrackingMode == kMcbm) {
        sizeX = 1.5;
        sizeY = 1.5;
        }
        hResidualsBeforeAlignmentX.push_back( hResidualsBeforeAlignmentX.push_back(
        new TH1F(Form("ResBeforeAli%s_X", n1), Form("X-Residuals Before Alignment%s", n2), 100, -1.5, 1.5)); new TH1F(Form("ResBeforeAli%s_X", n1), Form("X-Residuals Before Alignment%s", n2), 100, -sizeX, sizeX));
        hResidualsBeforeAlignmentY.push_back( hResidualsBeforeAlignmentY.push_back(
        new TH1F(Form("ResBeforeAli%s_Y", n1), Form("Y-Residuals Before Alignment%s", n2), 100, -1.5, 1.5)); new TH1F(Form("ResBeforeAli%s_Y", n1), Form("Y-Residuals Before Alignment%s", n2), 100, -sizeY, sizeY));
        hResidualsAfterAlignmentX.push_back( hResidualsAfterAlignmentX.push_back(
        new TH1F(Form("ResAfterAli%s_X", n1), Form("X-Residuals After Alignment%s", n2), 100, -1.5, 1.5)); new TH1F(Form("ResAfterAli%s_X", n1), Form("X-Residuals After Alignment%s", n2), 100, -sizeX, sizeX));
        hResidualsAfterAlignmentY.push_back( hResidualsAfterAlignmentY.push_back(
        new TH1F(Form("ResAfterAli%s_Y", n1), Form("Y-Residuals After Alignment%s", n2), 100, -1.5, 1.5)); new TH1F(Form("ResAfterAli%s_Y", n1), Form("Y-Residuals After Alignment%s", n2), 100, -sizeY, sizeY));
        hPullsBeforeAlignmentX.push_back( hPullsBeforeAlignmentX.push_back(
        new TH1F(Form("PullBeforeAli%s_X", n1), Form("X-Pulls Before Alignment%s", n2), 100, -10., 10.)); new TH1F(Form("PullBeforeAli%s_X", n1), Form("X-Pulls Before Alignment%s", n2), 100, -10., 10.));
        ...@@ -208,11 +234,9 @@ void CbmBbaAlignmentTask::Exec(Option_t* /*opt*/) ...@@ -208,11 +234,9 @@ void CbmBbaAlignmentTask::Exec(Option_t* /*opt*/)
        ca::Framework* l1 = CbmL1::Instance()->fpAlgo; ca::Framework* l1 = CbmL1::Instance()->fpAlgo;
        const ca::Parameters<ca::fvec>& l1Par = l1->GetParameters(); const ca::Parameters<ca::fvec>& l1Par = l1->GetParameters();
        static int statGlobalTracks = 0;
        statGlobalTracks += fInputGlobalTracks->GetEntriesFast();
        // select tracks for alignment and store them // select tracks for alignment and store them
        if (fTrackingMode == kSts && fInputStsTracks) { if (fTrackingMode == kSts && fInputStsTracks) {
        for (int iTr = 0; iTr < fInputStsTracks->GetEntriesFast(); iTr++) { for (int iTr = 0; iTr < fInputStsTracks->GetEntriesFast(); iTr++) {
        if (static_cast<int>(fTracks.size()) >= fMaxNtracks) { if (static_cast<int>(fTracks.size()) >= fMaxNtracks) {
        break; break;
        ...@@ -222,12 +246,24 @@ void CbmBbaAlignmentTask::Exec(Option_t* /*opt*/) ...@@ -222,12 +246,24 @@ void CbmBbaAlignmentTask::Exec(Option_t* /*opt*/)
        TrackContainer t; TrackContainer t;
        if (!fFitter.CreateMvdStsTrack(t.fUnalignedTrack, iTr)) continue; if (!fFitter.CreateMvdStsTrack(t.fUnalignedTrack, iTr)) continue;
        t.MakeConsistent(); t.MakeConsistent();
        t.fUnalignedTrack.fMsQp0 = 1. / 0.1;
        t.fUnalignedTrack.fIsMsQp0Set = false;
        for (auto& n : t.fUnalignedTrack.fNodes) {
        n.fIsTimeSet = false;
        }
        fFitter.FitTrack(t.fUnalignedTrack); 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; const auto& par = t.fUnalignedTrack.fNodes[t.fUnalignedTrack.fFirstHitNode].fTrack;
        t.fUnalignedTrack.fMsQp0 = par.GetQp()[0];
        //t.fUnalignedTrack.fIsMsQp0Set = true;
        if (t.fNstsHits < l1Par.GetNstationsActive(ca::EDetectorID::kSts)) continue;
        if (!std::isfinite((ca::fscal) par.GetQp()[0])) continue; if (!std::isfinite((ca::fscal) par.GetQp()[0])) continue;
        if (fabs(par.GetQp()[0]) > 1.) continue; // select tracks with min 1 GeV momentum if (fabs(par.GetQp()[0]) > 1.) continue; // select tracks with min 1 GeV momentum
        t.fAlignedTrack = t.fUnalignedTrack;
        fTracks.push_back(t); fTracks.push_back(t);
        } }
        } }
        ...@@ -251,6 +287,10 @@ void CbmBbaAlignmentTask::Exec(Option_t* /*opt*/) ...@@ -251,6 +287,10 @@ void CbmBbaAlignmentTask::Exec(Option_t* /*opt*/)
        t.fUnalignedTrack.fIsMsQp0Set = true; t.fUnalignedTrack.fIsMsQp0Set = true;
        t.fUnalignedTrack.fNodes[t.fUnalignedTrack.fFirstHitNode].fTrack.Qp() = 0.; t.fUnalignedTrack.fNodes[t.fUnalignedTrack.fFirstHitNode].fTrack.Qp() = 0.;
        for (auto& n : t.fUnalignedTrack.fNodes) {
        n.fIsTimeSet = false;
        }
        fFitter.FitTrack(t.fUnalignedTrack); fFitter.FitTrack(t.fUnalignedTrack);
        t.fAlignedTrack = t.fUnalignedTrack; t.fAlignedTrack = t.fUnalignedTrack;
        //if (t.fNstsHits < 1) continue; //if (t.fNstsHits < 1) continue;
        ...@@ -260,8 +300,18 @@ void CbmBbaAlignmentTask::Exec(Option_t* /*opt*/) ...@@ -260,8 +300,18 @@ void CbmBbaAlignmentTask::Exec(Option_t* /*opt*/)
        } }
        } }
        LOG(info) << "Bba: " << fInputGlobalTracks->GetEntriesFast() << " global tracks in event, " << statGlobalTracks TClonesArray* inputTracks = nullptr;
        << " global tracks in total, " << fTracks.size() << " tracks selected for alignment"; if (fTrackingMode == kSts) {
        inputTracks = fInputStsTracks;
        }
        else if (fTrackingMode == kMcbm) {
        inputTracks = fInputGlobalTracks;
        }
        static int statTracks = 0;
        statTracks += inputTracks->GetEntriesFast();
        LOG(info) << "Bba: " << inputTracks->GetEntriesFast() << " tracks in event, " << statTracks << " tracks in total, "
        << fTracks.size() << " tracks selected for alignment";
        } }
        ...@@ -300,20 +350,44 @@ double CbmBbaAlignmentTask::CostFunction(const std::vector<double>& par) ...@@ -300,20 +350,44 @@ double CbmBbaAlignmentTask::CostFunction(const std::vector<double>& par)
        ApplyAlignment(par); ApplyAlignment(par);
        fChi2Total = 0.; std::vector<double> chi2Thread(fNthreads, 0.);
        fNdfTotal = 0; std::vector<long> ndfThread(fNthreads, 0);
        for (auto& t : fTracks) { auto fitThread = [&](int iThread) {
        if (!t.fIsActive) continue; CbmKfTrackFitter fit(fFitter);
        fFitter.FitTrack(t.fAlignedTrack); for (unsigned int itr = iThread; itr < fTracks.size(); itr += fNthreads) {
        double chi2 = t.fAlignedTrack.fNodes[t.fAlignedTrack.fFirstHitNode].fTrack.GetChiSq()[0]; auto& t = fTracks[itr];
        double ndf = t.fAlignedTrack.fNodes[t.fAlignedTrack.fFirstHitNode].fTrack.GetNdf()[0]; if (!t.fIsActive) continue;
        if (!std::isfinite(chi2) || chi2 <= 0. || ndf <= 0.) { fit.FitTrack(t.fAlignedTrack);
        chi2 = 0.; double chi2 = t.fAlignedTrack.fNodes[t.fAlignedTrack.fFirstHitNode].fTrack.GetChiSq()[0];
        ndf = 0.; double ndf = t.fAlignedTrack.fNodes[t.fAlignedTrack.fFirstHitNode].fTrack.GetNdf()[0];
        if (!std::isfinite(chi2) || chi2 <= 0. || ndf <= 0.) {
        chi2 = 0.;
        ndf = 0.;
        }
        chi2Thread[iThread] += chi2;
        ndfThread[iThread] += (int) ndf;
        } }
        fChi2Total += chi2; };
        fNdfTotal += (int) ndf;
        std::vector<std::thread> threads(fNthreads);
        // run n threads
        for (int i = 0; i < fNthreads; i++) {
        threads[i] = std::thread(fitThread, i);
        }
        // wait for the threads to finish
        for (auto& th : threads) {
        th.join();
        }
        fChi2Total = 0.;
        fNdfTotal = 0;
        for (int i = 0; i < fNthreads; i++) {
        fChi2Total += chi2Thread[i];
        fNdfTotal += ndfThread[i];
        } }
        double cost = fChi2Total / (fNdfTotal + 1); double cost = fChi2Total / (fNdfTotal + 1);
        ...@@ -351,34 +425,26 @@ void CbmBbaAlignmentTask::Finish() ...@@ -351,34 +425,26 @@ void CbmBbaAlignmentTask::Finish()
        bba::Parameter& p = par[ip]; bba::Parameter& p = par[ip];
        p.SetActive(0); p.SetActive(0);
        p.SetValue(0.); p.SetValue(0.);
        p.SetBoundaryMin(-4.); // +-1 cm alignment range p.SetBoundaryMin(-20.); // +-1 cm alignment range
        p.SetBoundaryMax(4.); p.SetBoundaryMax(20.);
        p.SetStepMin(1.e-4); p.SetStepMin(1.e-6);
        p.SetStepMax(0.1); p.SetStepMax(0.1);
        p.SetStep(1.e-2); p.SetStep(1.e-2);
        } }
        par[0].SetActive(0); // fix the first station // set active parameters
        par[1].SetActive(0);
        par[2].SetActive(0);
        par[3 * (fNalignmentBodies - 1) + 0].SetActive(0); // fix the last station
        par[3 * (fNalignmentBodies - 1) + 1].SetActive(0);
        par[3 * (fNalignmentBodies - 1) + 2].SetActive(0);
        // set active parameters for other stations
        for (int i = 1; i < fNalignmentBodies - 1; i++) {
        par[3 * i + 0].SetActive(0);
        par[3 * i + 1].SetActive(0);
        par[3 * i + 2].SetActive(0);
        }
        for (int i = 1; i < fNalignmentBodies - 1; i++) { for (int i = 1; i < fNalignmentBodies - 1; i++) {
        //for (int i = 1; i < 2; i++) { if (fTrackingMode == kSts && i == fNalignmentBodies / 2) {
        par[3 * i + 0].SetActive(1); par[3 * i + 0].SetActive(0);
        par[3 * i + 1].SetActive(1); par[3 * i + 1].SetActive(1);
        par[3 * i + 2].SetActive(1); par[3 * i + 2].SetActive(0);
        }
        else {
        par[3 * i + 0].SetActive(1);
        par[3 * i + 1].SetActive(1);
        par[3 * i + 2].SetActive(1);
        }
        } }
        gRandom->SetSeed(1); gRandom->SetSeed(1);
        ...@@ -388,8 +454,8 @@ void CbmBbaAlignmentTask::Finish() ...@@ -388,8 +454,8 @@ void CbmBbaAlignmentTask::Finish()
        bba::Parameter& py = par[3 * is + 1]; bba::Parameter& py = par[3 * is + 1];
        bba::Parameter& pz = par[3 * is + 2]; bba::Parameter& pz = par[3 * is + 2];
        // +- 0.5 cm misalignment // +- d cm misalignment
        double d = .2; double d = 1.;
        if (px.IsActive()) { if (px.IsActive()) {
        px.SetValue(gRandom->Uniform(2. * d) - d); px.SetValue(gRandom->Uniform(2. * d) - d);
        } }
        ...@@ -432,13 +498,12 @@ void CbmBbaAlignmentTask::Finish() ...@@ -432,13 +498,12 @@ void CbmBbaAlignmentTask::Finish()
        t.fIsActive = false; t.fIsActive = false;
        } }
        } }
        std::cout << "recalculate initial cost function..." << std::endl; std::cout << "reject bad tracks and recalculate the initial cost function..." << std::endl;
        fCostInitial = CostFunction(parMisaligned); fCostInitial = CostFunction(parMisaligned);
        fFixedNdf = -1; //fNdfTotal; fFixedNdf = -1; //fNdfTotal;
        for (auto& t : fTracks) { for (auto& t : fTracks) {
        std::cout << "THIS IS THE HISTOGRAM LOOP" << std::endl;
        if (!t.fIsActive) continue; if (!t.fIsActive) continue;
        // calculate the residuals for all tracks at each fitNode before alignment // calculate the residuals for all tracks at each fitNode before alignment
        for (unsigned int in = 0; in < t.fAlignedTrack.fNodes.size(); in++) { for (unsigned int in = 0; in < t.fAlignedTrack.fNodes.size(); in++) {
        ...@@ -464,31 +529,40 @@ void CbmBbaAlignmentTask::Finish() ...@@ -464,31 +529,40 @@ void CbmBbaAlignmentTask::Finish()
        } }
        } }
        std::cout << "ideal cost function..." << std::endl; std::cout << "calculate ideal cost function..." << std::endl;
        fCostIdeal = CostFunction(parIdeal); fCostIdeal = CostFunction(parIdeal);
        LOG(info) << "Initial nTracks: " << fTracks.size() << " cost " << fCostInitial; LOG(info) << " Cost function for the true parameters: " << fCostIdeal;
        LOG(info) << " cost function for the true parameters: " << fCostIdeal; LOG(info) << " Cost function for the initial parameters: " << fCostInitial;
        alignment.align(); alignment.align();
        LOG(info) << " cost function for the true parameters: " << fCostIdeal; double costResult = CostFunction(alignment.getResult());
        LOG(info) << " Misaligned parameters: "; CostFunction(alignment.getResult());
        CostFunction(alignment.getResult());
        LOG(info) << " Cost function for the true parameters: " << fCostIdeal;
        LOG(info) << " Cost function for the initial parameters: " << fCostInitial;
        LOG(info) << " Cost function for the aligned parameters: " << costResult;
        LOG(info) << "\n Misaligned parameters: ";
        for (int is = 0; is < fNalignmentBodies; is++) { for (int is = 0; is < fNalignmentBodies; is++) {
        const std::vector<double>& r = parMisaligned; 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) << "Module " << is << ": x " << r[3 * is + 0] << " y " << r[3 * is + 1] << " z " << r[3 * is + 2];
        } }
        LOG(info) << " Alignment results: "; LOG(info) << "\n Alignment results: ";
        for (int is = 0; is < fNalignmentBodies; is++) { for (int is = 0; is < fNalignmentBodies; is++) {
        const std::vector<double>& r = alignment.getResult(); 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]; LOG(info) << "Module " << is << ": x " << r[3 * is + 0] << " y " << r[3 * is + 1] << " z " << r[3 * is + 2];
        } }
        CostFunction(alignment.getResult()); LOG(info) << "\n";
        LOG(info) << "Bba: " << fTracks.size() << " tracks used for alignment";
        for (auto& t : fTracks) { for (auto& t : fTracks) {
        if (!t.fIsActive) continue; if (!t.fIsActive) continue;
        ...@@ -505,10 +579,13 @@ void CbmBbaAlignmentTask::Finish() ...@@ -505,10 +579,13 @@ void CbmBbaAlignmentTask::Finish()
        hResidualsAfterAlignmentX[node.fReference1 + 1]->Fill(x_residual); hResidualsAfterAlignmentX[node.fReference1 + 1]->Fill(x_residual);
        hResidualsAfterAlignmentY[node.fReference1 + 1]->Fill(y_residual); hResidualsAfterAlignmentY[node.fReference1 + 1]->Fill(y_residual);
        hPullsAfterAlignmentX[0]->Fill(x_residual); double x_pull = x_residual / sqrt(node.fMxy.Dx2()[0] + node.fTrack.GetCovariance(0, 0)[0]);
        hPullsAfterAlignmentY[0]->Fill(y_residual); double y_pull = y_residual / sqrt(node.fMxy.Dy2()[0] + node.fTrack.GetCovariance(1, 1)[0]);
        hPullsAfterAlignmentX[node.fReference1 + 1]->Fill(x_residual);
        hPullsAfterAlignmentY[node.fReference1 + 1]->Fill(y_residual); hPullsAfterAlignmentX[0]->Fill(x_pull);
        hPullsAfterAlignmentY[0]->Fill(y_pull);
        hPullsAfterAlignmentX[node.fReference1 + 1]->Fill(x_pull);
        hPullsAfterAlignmentY[node.fReference1 + 1]->Fill(y_pull);
        } }
        } }
        ......
        reco/alignment/CbmBbaAlignmentTask.h
        + 2
        − 0
        View file @ 2ba38019
        ...@@ -96,6 +96,8 @@ class CbmBbaAlignmentTask : public FairTask { ...@@ -96,6 +96,8 @@ class CbmBbaAlignmentTask : public FairTask {
        TClonesArray* fInputGlobalTrackMatches{nullptr}; // Mc info for debugging TClonesArray* fInputGlobalTrackMatches{nullptr}; // Mc info for debugging
        TClonesArray* fInputStsTrackMatches{nullptr}; // Mc info for debugging TClonesArray* fInputStsTrackMatches{nullptr}; // Mc info for debugging
        int fNthreads = 1;
        CbmKfTrackFitter fFitter; CbmKfTrackFitter fFitter;
        // collection of selected tracks and hits // collection of selected tracks and hits
        ......
        0% or .
        You are about to add 0 people to the discussion. Proceed with caution.
        Finish editing this message first!
        Please register or to comment