: mDetMask(detmask), mMPsrc(trcmask), mUseMC(useMC), mInstanceID(instID)

{

setCosmic(cosmic);

init();

{

if (mDetMask[DetID::ITS]) {

addDetector(new AlignableDetectorITS(this));

}

if (mDetMask[DetID::TRD]) {

addDetector(new AlignableDetectorTRD(this));

}

if (mDetMask[DetID::TPC]) {

addDetector(new AlignableDetectorTPC(this));

}

if (mDetMask[DetID::TOF]) {

addDetector(new AlignableDetectorTOF(this));

}

for (int src = GIndex::NSources; src--;) {

if (mMPsrc[src]) {

mTrackSources.push_back(src);

}

}

mVtxSens = std::make_unique<EventVertex>(this);

mVtxSens->setInternalID(1);

const auto& algConf = AlignConfig::Instance();

if (algConf.MPRecOutFraction > 0. || mInstanceID == 0) {

mMPRecOutFraction = std::abs(algConf.MPRecOutFraction);

}

if (algConf.controlFraction > 0. || mInstanceID == 0) {

mControlFraction = std::abs(algConf.controlFraction);

}

{

o2::steer::MCKinematicsReader mcReader;

if (mUseMC) {

if (!mcReader.initFromDigitContext("collisioncontext.root")) {

throw std::invalid_argument("initialization of MCKinematicsReader failed");

}

}

auto timerStart = std::chrono::system_clock::now();

int nVtx = 0, nVtxAcc = 0, nTrc = 0, nTrcAcc = 0;

for (auto id = DetID::First; id <= DetID::Last; id++) {

auto* det = getDetector(id);

if (det) {

det->prepareDetectorData(); // in case the detector needs to preprocess the RecoContainer data

}

}

auto primVertices = mRecoData->getPrimaryVertices();

auto primVer2TRefs = mRecoData->getPrimaryVertexMatchedTrackRefs();

auto primVerGIs = mRecoData->getPrimaryVertexMatchedTracks();

const auto& algConf = AlignConfig::Instance();

// process vertices with contributor tracks

std::unordered_map<GIndex, bool> ambigTable;

int nvRefs = primVer2TRefs.size();

bool fieldON = std::abs(PropagatorD::Instance()->getNominalBz()) > 0.1;

for (int ivref = 0; ivref < nvRefs; ivref++) {

const o2::dataformats::PrimaryVertex* vtx = (ivref < nvRefs - 1) ? &primVertices[ivref] : nullptr;

bool useVertexConstrain = false;

if (vtx) {

auto nContrib = vtx->getNContributors();

// check cov matrix since data reconstructed with < 6797a257f5ab8ffaec32d56dddb0a321939bdf1c may have negative errors

if (vtx->getSigmaX2() < 0. || vtx->getSigmaY2() < 0. || vtx->getSigmaZ2() < 0.) {

continue;

}

useVertexConstrain = nContrib >= algConf.vtxMinCont && nContrib <= algConf.vtxMaxCont;

mStat.data[ProcStat::kInput][ProcStat::kVertices]++;

}

auto& trackRef = primVer2TRefs[ivref];

if (algConf.verbose > 1) {

LOGP(info, "processing vtref {} of {} with {} tracks, {}", ivref, nvRefs, trackRef.getEntries(), vtx ? vtx->asString() : std::string{});

}

nVtx++;

bool newVtx = true;

for (int src : mTrackSources) {

if ((GIndex::getSourceDetectorsMask(src) & mDetMask).none()) { // do we need this source?

continue;

}

int start = trackRef.getFirstEntryOfSource(src), end = start + trackRef.getEntriesOfSource(src);

for (int ti = start; ti < end; ti++) {

auto trackIndex = primVerGIs[ti];

mAlgTrack->setCurrentTrackID(trackIndex);

bool tpcIn = false;

if (trackIndex.isAmbiguous()) {

auto& ambSeen = ambigTable[trackIndex];

if (ambSeen) { // processed

continue;

}

ambSeen = true;

}

mStat.data[ProcStat::kInput][ProcStat::kTracks]++;

if (vtx) {

mStat.data[ProcStat::kInput][ProcStat::kTracksWithVertex]++;

}

int npnt = 0;

auto contributorsGID = mRecoData->getSingleDetectorRefs(trackIndex);

std::string trComb;

for (int ig = 0; ig < GIndex::NSources; ig++) {

if (contributorsGID[ig].isIndexSet()) {

trComb += " " + contributorsGID[ig].asString();

}

}

if (algConf.verbose > 1) {

LOG(info) << "processing track " << trackIndex.asString() << " contributors: " << trComb;

}

resetForNextTrack();

nTrc++;

// RS const auto& trcOut = mRecoData->getTrackParamOut(trackIndex);

auto trcOut = mRecoData->getTrackParamOut(trackIndex);

const auto& trcIn = mRecoData->getTrackParam(trackIndex);

// check detectors contributions

AlignableDetector* det = nullptr;

int ndet = 0, npntDet = 0;

if ((det = getDetector(DetID::ITS))) {

if (contributorsGID[GIndex::ITS].isIndexSet() && (npntDet = det->processPoints(contributorsGID[GIndex::ITS], algConf.minITSClusters, false)) > 0) {

npnt += npntDet;

ndet++;

} else if (mAllowAfterburnerTracks && contributorsGID[GIndex::ITSAB].isIndexSet() && (npntDet = det->processPoints(contributorsGID[GIndex::ITSAB], 2, false)) > 0) {

npnt += npntDet;

ndet++;

} else {

continue;

}

}

if ((det = getDetector(DetID::TPC)) && contributorsGID[GIndex::TPC].isIndexSet()) {

float t0 = 0, t0err = 0;

mRecoData->getTrackTime(trackIndex, t0, t0err);

((AlignableDetectorTPC*)det)->setTrackTimeStamp(t0);

npntDet = det->processPoints(contributorsGID[GIndex::TPC], algConf.minTPCClusters, false);

if (npntDet > 0) {

npnt += npntDet;

ndet++;

tpcIn = true;

}

}

if ((det = getDetector(DetID::TRD)) && contributorsGID[GIndex::TRD].isIndexSet() && (npntDet = det->processPoints(contributorsGID[GIndex::TRD], algConf.minTRDTracklets, false)) > 0) {

npnt += npntDet;

ndet++;

}

if ((det = getDetector(DetID::TOF)) && contributorsGID[GIndex::TOF].isIndexSet() && (npntDet = det->processPoints(contributorsGID[GIndex::TOF], algConf.minTOFClusters, false)) > 0) {

npnt += npntDet;

ndet++;

}

// other detectors

if (algConf.verbose > 1) {

LOGP(info, "processing track {} {} of vtref {}, Ndets:{}, Npoints: {}, use vertex: {} | Kin: {} Kout: {}", ti, trackIndex.asString(), ivref, ndet, npnt, useVertexConstrain && trackIndex.isPVContributor(), trcIn.asString(), trcOut.asString());

}

if (ndet < algConf.minDetectors || (tpcIn && ndet == 1)) { // we don't want TPC only track

continue;

}

if (npnt < algConf.minPointTotal) {

if (algConf.verbose > 0) {

LOGP(info, "too few points {} < {}", npnt, algConf.minPointTotal);

}

continue;

}

bool vtxCont = false;

if (trackIndex.isPVContributor() && useVertexConstrain) {

mAlgTrack->copyFrom(trcIn); // copy kinematices of inner track just for propagation to the vertex

if (addVertexConstraint(*vtx)) {

mAlgTrack->setRefPoint(mRefPoint.get()); // set vertex as a reference point

vtxCont = true;

}

}

mAlgTrack->copyFrom(trcOut); // copy kinematices of outer track as the refit will be done inward

mAlgTrack->setFieldON(fieldON);

mAlgTrack->sortPoints();

int pntMeas = mAlgTrack->getInnerPointID() - 1;

if (pntMeas < 0) { // this should not happen

mAlgTrack->Print("p meas");

LOG(error) << "AliAlgTrack->GetInnerPointID() cannot be 0";

}

if (!mAlgTrack->iniFit()) {

if (algConf.verbose > 0) {

LOGP(warn, "iniFit failed");

}

continue;

}

// compare refitted and original track

if (mDebugOutputLevel) {

trackParam_t trcAlgRef(*mAlgTrack.get());

std::array<double, 5> dpar{};

std::array<double, 15> dcov{};

for (int i = 0; i < 5; i++) {

dpar[i] = trcIn.getParam(i);

}

for (int i = 0; i < 15; i++) {

dcov[i] = trcIn.getCov()[i];

}

trackParam_t trcOrig(trcIn.getX(), trcIn.getAlpha(), dpar, dcov, trcIn.getCharge());

if (PropagatorD::Instance()->propagateToAlphaX(trcOrig, trcAlgRef.getAlpha(), trcAlgRef.getX(), true)) {

(*mDBGOut) << "trcomp"

<< "orig=" << trcOrig << "fit=" << trcAlgRef << "\n";

}

}

// RS: this is to substitute the refitter track by MC truth, just for debugging

/*

if (!mAlgTrack->processMaterials()) {

if (algConf.verbose > 0) {

LOGP(warn, "processMaterials failed");

}

continue;

}

mAlgTrack->defineDOFs();

if (!mAlgTrack->calcResidDeriv()) {

if (algConf.verbose > 0) {

LOGP(warn, "calcResidDeriv failed");

}

continue;

}

if (mDebugOutputLevel && mAlgTrackDbg.setTrackParam(mAlgTrack.get())) {

mAlgTrackDbg.mGID = trackIndex;

(*mDBGOut) << "algtrack"

<< "runNumber=" << mTimingInfo.runNumber

<< "tfID=" << mTimingInfo.tfCounter

<< "orbit=" << mTimingInfo.firstTForbit

<< "bz=" << PropagatorD::Instance()->getNominalBz()

<< "t=" << mAlgTrackDbg << "\n";

}

if (mUseMC && mDebugOutputLevel > 1) {

auto lbl = mRecoData->getTrackMCLabel(trackIndex);

if (lbl.isValid()) {

std::vector<float> pntX, pntY, pntZ, trcX, trcY, trcZ, prpX, prpY, prpZ, alpha, xsens, pntXTF, pntYTF, pntZTF, resY, resZ;

std::vector<int> detid, volid;

o2::MCTrack mcTrack = *mcReader.getTrack(lbl);

trackParam_t recTrack{*mAlgTrack};

for (int ip = 0; ip < mAlgTrack->getNPoints(); ip++) {

double tmp[3], tmpg[3];

auto* pnt = mAlgTrack->getPoint(ip);

auto* sens = pnt->getSensor();

detid.emplace_back(pnt->getDetID());

volid.emplace_back(pnt->getVolID());

TGeoHMatrix t2g;

sens->getMatrixT2G(t2g);

t2g.LocalToMaster(pnt->getXYZTracking(), tmpg);

pntX.emplace_back(tmpg[0]);

pntY.emplace_back(tmpg[1]);

pntZ.emplace_back(tmpg[2]);

double xyz[3]{pnt->getXTracking(), pnt->getYTracking(), pnt->getZTracking()};

xyz[1] += mAlgTrack->getResidual(0, ip);

xyz[2] += mAlgTrack->getResidual(1, ip);

t2g.LocalToMaster(xyz, tmpg);

trcX.emplace_back(tmpg[0]);

trcY.emplace_back(tmpg[1]);

trcZ.emplace_back(tmpg[2]);

pntXTF.emplace_back(pnt->getXTracking());

pntYTF.emplace_back(pnt->getYTracking());

pntZTF.emplace_back(pnt->getZTracking());

resY.emplace_back(mAlgTrack->getResidual(0, ip));

resZ.emplace_back(mAlgTrack->getResidual(1, ip));

alpha.emplace_back(pnt->getAlphaSens());

xsens.emplace_back(pnt->getXSens());

}

(*mDBGOut) << "mccomp"

<< "mcTr=" << mcTrack << "recTr=" << recTrack << "gid=" << trackIndex << "lbl=" << lbl << "vtxConst=" << vtxCont

<< "pntX=" << pntX << "pntY=" << pntY << "pntZ=" << pntZ

<< "trcX=" << trcX << "trcY=" << trcY << "trcZ=" << trcZ

<< "alp=" << alpha << "xsens=" << xsens

<< "pntXTF=" << pntXTF << "pntYTF=" << pntYTF << "pntZTF=" << pntZTF

<< "resY=" << resY << "resZ=" << resZ

<< "detid=" << detid << "volid=" << volid << "\n";

}

}

mStat.data[ProcStat::kAccepted][ProcStat::kTracks]++;

if (vtxCont) {

mStat.data[ProcStat::kAccepted][ProcStat::kTracksWithVertex]++;

}

nTrcAcc++;

if (newVtx) {

newVtx = false;

mStat.data[ProcStat::kAccepted][ProcStat::kVertices]++;

nVtxAcc++;

}

storeProcessedTrack(trackIndex);

}

}

}

auto timerEnd = std::chrono::system_clock::now();

std::chrono::duration<float, std::milli> duration = timerEnd - timerStart;

LOGP(info, "Processed TF {}: {} vertices ({} used), {} tracks ({} used) in {} ms", mNTF, nVtx, nVtxAcc, nTrc, nTrcAcc, duration.count());

mNTF++;

{

auto timerStart = std::chrono::system_clock::now();

const auto tracks = mRecoData->getCosmicTracks();

if (!tracks.size()) {

LOGP(info, "Skipping TF {}: No cosmic tracks", mNTF);

mNTF++;

return;

}

int nTrc = 0, nTrcAcc = 0;

for (auto id = DetID::First; id <= DetID::Last; id++) {

auto* det = getDetector(id);

if (det) {

det->prepareDetectorData(); // in case the detector needs to preprocess the RecoContainer data

}

}

const auto& algConf = AlignConfig::Instance();

bool fieldON = std::abs(PropagatorD::Instance()->getNominalBz()) > 0.1;

for (const auto& track : tracks) {

resetForNextTrack();

nTrc++;

mStat.data[ProcStat::kInput][ProcStat::kCosmic]++;

std::array<GTrackID, GTrackID::NSources> contributorsGID[2] = {mRecoData->getSingleDetectorRefs(track.getRefBottom()), mRecoData->getSingleDetectorRefs(track.getRefTop())};

bool hasTRD = false, hasITS = false, hasTPC = false, hasTOF = false;

if (contributorsGID[0][GTrackID::TRD].isIndexSet() || contributorsGID[1][GTrackID::TRD].isIndexSet()) {

hasTRD = true;

}

if (contributorsGID[0][GTrackID::TOF].isIndexSet() || contributorsGID[1][GTrackID::TOF].isIndexSet()) {

hasTOF = true;

}

if (contributorsGID[0][GTrackID::TPC].isIndexSet() || contributorsGID[1][GTrackID::TPC].isIndexSet()) {

hasTPC = true;

}

if (contributorsGID[0][GTrackID::ITS].isIndexSet() || contributorsGID[1][GTrackID::ITS].isIndexSet()) {

hasITS = true;

}

// check detectors contributions

AlignableDetector* det = nullptr;

int ndet = 0, npnt = 0;

bool accTrack = true;

bool tpcIn = false;

if ((det = getDetector(DetID::ITS))) {

int npntDet = 0;

for (int ibt = 0; ibt < 2; ibt++) {

int npntDetBT = 0;

mAlgTrack->setCurrentTrackID(ibt ? track.getRefBottom() : track.getRefTop());

if (contributorsGID[ibt][GIndex::ITS].isIndexSet() && (npntDetBT = det->processPoints(contributorsGID[ibt][GIndex::ITS], algConf.minITSClustersCosmLeg, ibt)) < 0) {

accTrack = false;

break;

}

npntDet += npntDetBT;

}

if (!accTrack || npntDet < algConf.minITSClustersCosm) {

continue;

}

if (npntDet) {

ndet++;

npnt += npntDet;

}

}

if ((det = getDetector(DetID::TPC))) {

int npntDet = 0;

((AlignableDetectorTPC*)det)->setTrackTimeStamp(track.getTimeMUS().getTimeStamp());

for (int ibt = 0; ibt < 2; ibt++) {

int npntDetBT = 0;

mAlgTrack->setCurrentTrackID(ibt ? track.getRefBottom() : track.getRefTop());

if (contributorsGID[ibt][GIndex::TPC].isIndexSet() && (npntDetBT = det->processPoints(contributorsGID[ibt][GIndex::TPC], algConf.minTPCClustersCosmLeg, ibt)) < 0) {

accTrack = false;

break;

}

npntDet += npntDetBT;

}

if (!accTrack || npntDet < algConf.minTPCClustersCosm) {

continue;

}

if (npntDet) {

ndet++;

npnt += npntDet;

tpcIn = true;

}

}

if ((det = getDetector(DetID::TRD))) {

int npntDet = 0;

for (int ibt = 0; ibt < 2; ibt++) {

int npntDetBT = 0;

mAlgTrack->setCurrentTrackID(ibt ? track.getRefBottom() : track.getRefTop());

if (contributorsGID[ibt][GIndex::TRD].isIndexSet() && (npntDetBT = det->processPoints(contributorsGID[ibt][GIndex::TRD], algConf.minTRDTrackletsCosmLeg, ibt)) < 0) {

accTrack = false;

break;

}

npntDet += npntDetBT;

}

if (!accTrack || npntDet < algConf.minTRDTrackletsCosm) {

continue;

}

if (npntDet) {

ndet++;

npnt += npntDet;

}

}

if ((det = getDetector(DetID::TOF))) {

int npntDet = 0;

for (int ibt = 0; ibt < 2; ibt++) {

int npntDetBT = 0;

mAlgTrack->setCurrentTrackID(ibt ? track.getRefBottom() : track.getRefTop());

if (contributorsGID[ibt][GIndex::TOF].isIndexSet() && (npntDetBT = det->processPoints(contributorsGID[ibt][GIndex::TOF], algConf.minTOFClustersCosmLeg, ibt)) < 0) {

accTrack = false;

break;

}

npntDet += npntDetBT;

}

if (!accTrack || npntDet < algConf.minTOFClustersCosm) {

continue;

}

if (npntDet) {

ndet++;

npnt += npntDet;

}

}

if (algConf.verbose > 1) {

LOGP(info, "processing cosmic track B-Leg:{} T-Leg:{}, Ndets:{}, Npoints: {}", track.getRefBottom().asString(), track.getRefTop().asString(), ndet, npnt);

}

if (ndet < algConf.minDetectorsCosm /* || (tpcIn && ndet == 1)*/) {

continue;

}

if (npnt < algConf.minPointTotalCosm) {

if (algConf.verbose > 0) {

LOGP(info, "too few points {} < {}", npnt, algConf.minPointTotalCosm);

}

continue;

}

mAlgTrack->setCosmic(true);

mAlgTrack->copyFrom(mRecoData->getTrackParamOut(track.getRefBottom())); // copy kinematices of outer track as the refit will be done inward

mAlgTrack->setFieldON(fieldON);

mAlgTrack->sortPoints();

if (!mAlgTrack->iniFit()) {

if (algConf.verbose > 0) {

LOGP(warn, "iniFit failed");

}

continue;

}

if (!mAlgTrack->processMaterials()) {

if (algConf.verbose > 0) {

LOGP(warn, "processMaterials failed");

}

continue;

}

mAlgTrack->defineDOFs();

if (!mAlgTrack->calcResidDeriv()) {

if (algConf.verbose > 0) {

LOGP(warn, "calcResidDeriv failed");

}

continue;

}

if (mDebugOutputLevel && mAlgTrackDbg.setTrackParam(mAlgTrack.get())) {

mAlgTrackDbg.mGID = track.getRefBottom();

mAlgTrackDbg.mGIDCosmUp = track.getRefTop();

(*mDBGOut) << "algtrack"

<< "runNumber=" << mTimingInfo.runNumber

<< "tfID=" << mTimingInfo.tfCounter

<< "orbit=" << mTimingInfo.firstTForbit

<< "bz=" << PropagatorD::Instance()->getNominalBz()

<< "t=" << mAlgTrackDbg << "\n";

}

storeProcessedTrack();

mStat.data[ProcStat::kAccepted][ProcStat::kCosmic]++;

nTrcAcc++;

}

auto timerEnd = std::chrono::system_clock::now();

std::chrono::duration<float, std::milli> duration = timerEnd - timerStart;

LOGP(info, "Processed cosmic TF {}: {} tracks ({} used) in {} ms", mNTF, nTrc, nTrcAcc, duration.count());

mNTF++;

{

// init all detectors geometry

//

if (getInitGeomDone()) {

return;

}

//

mAlgTrack = std::make_unique<AlignmentTrack>();

mRefPoint = std::make_unique<AlignmentPoint>();

//

int dofCnt = 0;

// special fake sensor for vertex constraint point

// it has special T2L matrix adjusted for each track, no need to init it here

mVtxSens->prepareMatrixL2G();

mVtxSens->prepareMatrixL2GIdeal();

dofCnt += mVtxSens->getNDOFs();

//

for (auto id = DetID::First; id <= DetID::Last; id++) {

auto* det = getDetector(id);

if (det) {

dofCnt += det->initGeom();

}

}

if (!dofCnt) {

LOG(fatal) << "No DOFs found";

}

//

//

for (auto id = DetID::First; id <= DetID::Last; id++) {

auto* det = getDetector(id);

if (!det || det->isDisabled()) {

continue;

}

det->cacheReferenceCCDB();

}

//

assignDOFs();

LOG(info) << "Booked " << dofCnt << " global parameters";

//

setInitGeomDone();

//

{

// scan all free global parameters, link detectors to array of params

//

if (getInitDOFsDone()) {

LOG(info) << "initDOFs was already done, just reassigning " << getNDOFs() << "DOFs arrays/labels";

assignDOFs();

return;

}

const auto& conf = AlignConfig::Instance();

mNDOFs = 0;

int ndfAct = 0;

assignDOFs();

int nact = 0;

mVtxSens->initDOFs();

for (auto id = DetID::First; id <= DetID::Last; id++) {

AlignableDetector* det = getDetector(id);

if (det && !det->isDisabled()) {

det->initDOFs();

nact++;

ndfAct += det->getNDOFs();

}

}

for (int i = 0; i < NTrackTypes; i++) {

if (nact < conf.minDetAcc[i]) {

LOG(fatal) << nact << " detectors are active, while " << conf.minDetAcc[i] << " in track are asked";

}

}

LOG(info) << mNDOFs << " global parameters " << mNDet << " detectors, " << ndfAct << " in " << nact << " active detectors";

addAutoConstraints();

setInitDOFsDone();

{

// add parameters/labels arrays to volumes. If the Controller is read from the file, this method need

// to be called (of initDOFs should be called)

//

int ndfOld = -1;

if (mNDOFs > 0) {

ndfOld = mNDOFs;

}

mNDOFs = 0;

//

// reserve

int ndofTOT = mVtxSens->getNDOFs();

for (auto id = DetID::First; id <= DetID::Last; id++) {

AlignableDetector* det = getDetector(id);

if (!det) {

continue;

}

ndofTOT += det->getNDOFs();

}

mGloParVal.clear();

mGloParErr.clear();

mGloParLab.clear();

mLbl2ID.clear();

mGloParVal.reserve(ndofTOT);

mGloParErr.reserve(ndofTOT);

mGloParLab.reserve(ndofTOT);

mVtxSens->assignDOFs();

for (auto id = DetID::First; id <= DetID::Last; id++) {

AlignableDetector* det = getDetector(id);

if (!det) {

continue;

}

det->assignDOFs();

}

LOG(info) << "Assigned parameters/labels arrays for " << mNDOFs << " DOFs";

if (ndfOld > 0 && ndfOld != mNDOFs) {

LOG(error) << "Recalculated NDOFs=" << mNDOFs << " not equal to saved NDOFs=" << ndfOld;

}

// build Lbl -> parID table

for (int i = 0; i < ndofTOT; i++) {

int& id = mLbl2ID[mGloParLab[i]];

if (id != 0) {

LOGP(fatal, "parameters {} and {} share the same label {}", id - 1, i, mGloParLab[i]);

}

id = i + 1;

}

//

{

// add detector constructed externally to alignment framework

mDetectors[det->getDetID()].reset(det);

mNDet++;

{

//validate detector pattern

return ((patt & mObligatoryDetPattern[mTracksType]) == mObligatoryDetPattern[mTracksType]) &&

patt.count() >= AlignConfig::Instance().minDetAcc[mTracksType];

{

//validate detectors pattern according to number of selected points

int ndOK = 0;

for (auto id = DetID::First; id <= DetID::Last; id++) {

AlignableDetector* det = getDetector(id);

if (!det || det->isDisabled(mTracksType)) {

continue;

}

if (npsel[id] < det->getNPointsSel(mTracksType)) {

if (det->isObligatory(mTracksType)) {

return false;

}

continue;

}

ndOK++;

}

return ndOK >= AlignConfig::Instance().minDetAcc[mTracksType];

{

// write alignment track

bool res = true;

const auto& conf = AlignConfig::Instance();

if (conf.MilleOut) {

res &= fillMilleData();

}

float rnd = gRandom->Rndm();

if (mMPRecOutFraction > rnd) {

res &= fillMPRecData(tid);

}

if ((mControlFraction > rnd) && mAlgTrack->testLocalSolution()) {

res &= fillControlData(tid);

}

//

if (!res) {

LOGP(error, "storeProcessedTrack failed");

}

return res;

{

// store MP2 data in Mille format

if (!mMille) {

const auto& conf = AlignConfig::Instance();

mMilleFileName = fmt::format("{}_{:08d}_{:010d}{}", AlignConfig::Instance().mpDatFileName, mTimingInfo.runNumber, mTimingInfo.tfCounter, conf.MilleOutBin ? sMPDataExt : sMPDataTxtExt);

mMille = std::make_unique<Mille>(mMilleFileName.c_str(), conf.MilleOutBin);

}

if (!mAlgTrack->getDerivDone()) {

LOG(error) << "Track derivatives are not yet evaluated";

return false;

}

int np = mAlgTrack->getNPoints(), nDGloTot = 0; // total number global derivatives stored

int nParETP = mAlgTrack->getNLocExtPar(); // numnber of local parameters for reference track param

int nVarLoc = mAlgTrack->getNLocPar(); // number of local degrees of freedom in the track

//

const int* gloParID = mAlgTrack->getGloParID(); // IDs of global DOFs this track depends on

for (int ip = 0; ip < np; ip++) {

AlignmentPoint* pnt = mAlgTrack->getPoint(ip);

if (pnt->containsMeasurement()) {

int gloOffs = pnt->getDGloOffs(); // 1st entry of global derivatives for this point

int nDGlo = pnt->getNGloDOFs(); // number of global derivatives (number of DOFs it depends on)

if (!pnt->isStatOK()) {

pnt->incrementStat();

}

int milleIBufferG[nDGlo];

float milleDBufferG[nDGlo];

float milleDBufferL[nVarLoc];

std::memset(milleIBufferG, 0, sizeof(int) * nDGlo);

std::memset(milleDBufferG, 0, sizeof(float) * nDGlo);

std::memset(milleDBufferL, 0, sizeof(float) * nVarLoc);

// local der. array cannot be 0-suppressed by Mille construction, need to reset all to 0

for (int idim = 0; idim < 2; idim++) { // 2 dimensional orthogonal measurement

const double* deriv = mAlgTrack->getDResDLoc(idim, ip); // array of Dresidual/Dparams_loc

// derivatives over reference track parameters

for (int j = 0; j < nParETP; j++) {

milleDBufferL[j] = isZeroAbs(deriv[j]) ? 0. : deriv[j];

}

//

// point may depend on material variables within these limits

for (int j = pnt->getMinLocVarID(); j < pnt->getMaxLocVarID(); j++) {

milleDBufferL[j] = isZeroAbs(deriv[j]) ? 0. : deriv[j];

}

// derivatives over global params: this array can be 0-suppressed, no need to reset

int nGlo = 0;

deriv = mAlgTrack->getDResDGlo(idim, gloOffs);

const int* gloIDP(gloParID + gloOffs);

for (int j = 0; j < nDGlo; j++) {

milleDBufferG[nGlo] = isZeroAbs(deriv[j]) ? 0. : deriv[j]; // value of derivative

milleIBufferG[nGlo++] = getGloParLab(gloIDP[j]); // global DOF ID + 1 (Millepede needs positive labels)

}

mMille->mille(nVarLoc, milleDBufferL, nGlo, milleDBufferG, milleIBufferG, mAlgTrack->getResidual(idim, ip), Sqrt(pnt->getErrDiag(idim)));

nDGloTot += nGlo;

}

}

if (pnt->containsMaterial()) { // material point can add 4 or 5 otrhogonal pseudo-measurements

int nmatpar = pnt->getNMatPar(); // residuals (correction expectation value)

// const float* expMatCorr = pnt->getMatCorrExp(); // expected corrections (diagonalized)

const float* expMatCov = pnt->getMatCorrCov(); // their diagonalized error matrix

int offs = pnt->getMaxLocVarID() - nmatpar; // start of material variables

// here all derivatives are 1 = dx/dx

float milleDBufferL[nVarLoc];

std::memset(milleDBufferL, 0, sizeof(float) * nVarLoc);

for (int j = 0; j < nmatpar; j++) { // mat. "measurements" don't depend on global params

int j1 = j + offs;

milleDBufferL[j1] = 1.0; // only 1 non-0 derivative

// mMille->mille(nVarLoc,milleDBufferL,0, nullptr, nullptr, expMatCorr[j], Sqrt(expMatCov[j]));

// expectation for MS effect is 0

mMille->mille(nVarLoc, milleDBufferL, 0, nullptr, nullptr, 0, Sqrt(expMatCov[j]));

milleDBufferL[j1] = 0.0; // reset buffer

}

} // material "measurement"

} // loop over points

//

if (!nDGloTot) {

LOG(info) << "Track does not depend on free global parameters, discard";

mMille->clear();

return false;

}

mMille->finalise(); // store the record

return true;

{

// store MP2 in MPRecord format

if (!mMPRecFile) {

initMPRecOutput();

}

mMPRecord.clear();

if (!mMPRecord.fillTrack(*mAlgTrack.get(), mGloParLab)) {

return false;

}

mMPRecord.setRun(mRunNumber);

mMPRecord.setFirstTFOrbit(mTimingInfo.firstTForbit);

mMPRecord.setTrackID(tid);

mMPRecTree->Fill();

return true;

{

// store control residuals

if (!mResidFile) {

initResidOutput();

}

int nps, np = mAlgTrack->getNPoints();

nps = (!mRefPoint->containsMeasurement()) ? np - 1 : np; // ref point is dummy?

if (nps < 0) {

return true;

}

mCResid.clear();

if (!mCResid.fillTrack(*mAlgTrack.get(), AlignConfig::Instance().KalmanResid)) {

return false;

}

mCResid.setRun(mRunNumber);

mCResid.setFirstTFOrbit(mTimingInfo.firstTForbit);

mCResid.setBz(o2::base::PropagatorD::Instance()->getNominalBz());

mCResid.setTrackID(tid);

// if (isCosmic()) {

// mCResid.setInvTrackID(tid);

// }

mResidTree->Fill();

return true;

{

mTimingInfo = ti;

LOGP(info, "TIMING {} {}", ti.runNumber, ti.creation);

if (ti.runNumber != mRunNumber) {

mRunNumber = ti.runNumber;

}

{

// print info

TString opts = opt;

opts.ToLower();

printf("%5d DOFs in %d detectors\n", mNDOFs, mNDet);

if (getMPAlignDone()) {

printf("ALIGNMENT FROM MILLEPEDE SOLUTION IS APPLIED\n");

}

//

for (auto id = DetID::First; id <= DetID::Last; id++) {

AlignableDetector* det = getDetector(id);

if (!det) {

continue;

}

det->Print(opt);

}

if (!opts.IsNull()) {

printf("\nSpecial sensor for Vertex Constraint\n");

mVtxSens->Print(opt);

}

//

if (mRefRunNumber >= 0) {

printf("(%d)", mRefRunNumber);

}

AlignConfig::Instance().printKeyValues();

//

if (opts.Contains("stat")) {

printStatistics();

}

{

// print processing stat

mStat.print();

{

// reset detectors for next track

mRefPoint->clear();

mAlgTrack->Clear();

for (auto id = DetID::First; id <= DetID::Last; id++) {

AlignableDetector* det = getDetector(id);

if (det) {

det->reset();

}

}

{

// test track local solution

int npnt = mAlgTrack->getNPoints(), nlocpar = mAlgTrack->getNLocPar();

double mat[nlocpar][nlocpar], rhs[nlocpar];

std::memset(mat, 0, sizeof(double) * nlocpar * nlocpar);

std::memset(rhs, 0, sizeof(double) * nlocpar);

for (int ip = npnt; ip--;) {

AlignmentPoint* pnt = mAlgTrack->getPoint(ip);

if (pnt->containsMeasurement()) {

for (int idim = 2; idim--;) { // each point has 2 position residuals

double resid = mAlgTrack->getResidual(idim, ip), sg2inv = 1. / pnt->getErrDiag(idim); // residual and its inv. error

auto deriv = mAlgTrack->getDResDLoc(idim, ip); // array of Dresidual/Dparams

for (int parI = 0; parI < nlocpar; parI++) {

rhs[parI] -= deriv[parI] * resid * sg2inv;

for (int parJ = parI; parJ < nlocpar; parJ++) {

mat[parI][parJ] += deriv[parI] * deriv[parJ] * sg2inv;

}

}

} // loop over 2 orthogonal measurements at the point

} // derivarives at measured points

// if the point contains material, consider its expected kinks, eloss as measurements

if (pnt->containsMaterial()) { // at least 4 parameters: 2 spatial + 2 angular kinks with 0 expectaction

int npm = pnt->getNMatPar();

// const float* expMatCorr = pnt->getMatCorrExp(); // expected correction (diagonalized) // RS??