diff --git a/ALICE3/Core/CMakeLists.txt b/ALICE3/Core/CMakeLists.txt index 6db27c117d8..788711ac277 100644 --- a/ALICE3/Core/CMakeLists.txt +++ b/ALICE3/Core/CMakeLists.txt @@ -15,6 +15,14 @@ o2physics_add_library(ALICE3Core PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore) o2physics_target_root_dictionary(ALICE3Core - HEADERS TOFResoALICE3.h - DelphesO2TrackSmearer.h - LINKDEF ALICE3CoreLinkDef.h) + HEADERS TOFResoALICE3.h + DelphesO2TrackSmearer.h + LINKDEF ALICE3CoreLinkDef.h) + +o2physics_add_library(FastTracker + SOURCES FastTracker.cxx + PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore) + +o2physics_target_root_dictionary(FastTracker + HEADERS FastTracker.h + LINKDEF FastTrackerLinkDef.h) diff --git a/ALICE3/Core/DetLayer.h b/ALICE3/Core/DetLayer.h new file mode 100644 index 00000000000..7fa6aa25a53 --- /dev/null +++ b/ALICE3/Core/DetLayer.h @@ -0,0 +1,49 @@ +// Copyright 2019-2020 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. + +/// +/// \file DetLayer.h +/// \author David Dobrigkeit Chinellato +/// \since 11/03/2021 +/// \brief Basic struct to hold information regarding a detector layer to be used in fast simulation +/// + +#ifndef ALICE3_CORE_DETLAYER_H_ +#define ALICE3_CORE_DETLAYER_H_ + +#include "TString.h" + +namespace o2::fastsim +{ + +struct DetLayer { + // TString for holding name + TString name; + + // position variables + float r; // radius in centimeters + float z; // z dimension in centimeters + + // material variables + float x0; // radiation length + float xrho; // density + + // resolution variables for active layers + float resRPhi; // RPhi resolution in centimeters + float resZ; // Z resolution in centimeters + + // efficiency + float eff; // detection efficiency +}; + +} // namespace o2::fastsim + +#endif // ALICE3_CORE_DETLAYER_H_ diff --git a/ALICE3/Core/FastTracker.cxx b/ALICE3/Core/FastTracker.cxx new file mode 100644 index 00000000000..25acbba4af5 --- /dev/null +++ b/ALICE3/Core/FastTracker.cxx @@ -0,0 +1,380 @@ +// Copyright 2019-2020 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. + +#include +#include "TMath.h" +#include "TMatrixD.h" +#include "TRandom.h" +#include "TMatrixDSymEigen.h" +#include "FastTracker.h" + +namespace o2 +{ +namespace fastsim +{ + +// +-~-<*>-~-+-~-<*>-~-+-~-<*>-~-+-~-<*>-~-+-~-<*>-~-+-~-<*>-~-+-~-<*>-~-+-~-<*>-~-+ + +FastTracker::FastTracker() +{ + // base constructor + magneticField = 5; // in kiloGauss + applyZacceptance = false; + covMatFactor = 0.99f; + verboseLevel = 0; + covMatOK = 0; + covMatNotOK = 0; +} + +void FastTracker::AddLayer(TString name, float r, float z, float x0, float xrho, float resRPhi, float resZ, float eff) +{ + DetLayer newLayer{name.Data(), r, z, x0, xrho, resRPhi, resZ, eff}; + layers.push_back(newLayer); +} + +void FastTracker::Print() +{ + // print out layer setup + LOG(info) << "+-~-<*>-~-+-~-<*>-~-+-~-<*>-~-+-~-<*>-~-+-~-<*>-~-+-~-<*>-~-+-~-<*>-~-+"; + LOG(info) << " Printing detector layout with " << layers.size() << " effective elements: "; + for (uint32_t il = 0; il < layers.size(); il++) { + LOG(info) << " Layer #" << il << "\t" << layers[il].name.Data() << "\tr = " << Form("%.2f", layers[il].r) << "cm\tz = " << layers[il].z << "\t" + << "x0 = " << layers[il].x0 << "\txrho = " << layers[il].xrho << "\tresRPhi = " << layers[il].resRPhi << "\tresZ = " << layers[il].resZ << "\teff = " << layers[il].eff; + } + LOG(info) << "+-~-<*>-~-+-~-<*>-~-+-~-<*>-~-+-~-<*>-~-+-~-<*>-~-+-~-<*>-~-+-~-<*>-~-+"; +} + +void FastTracker::AddSiliconALICE3v4() +{ + LOG(info) << " Adding ALICE 3 v4 ITS layers"; + float x0IT = 0.001; // 0.1% + float x0OT = 0.005; // 0.5% + float xrhoIB = 1.1646e-02; // 50 mum Si + float xrhoOB = 1.1646e-01; // 500 mum Si + + float resRPhiIT = 0.00025; // 2.5 mum + float resZIT = 0.00025; // 2.5 mum + float resRPhiOT = 0.0005; // 5 mum + float resZOT = 0.0005; // 5 mum + float eff = 1.00; + + layers.push_back(DetLayer{"bpipe0", 0.48, 250, 0.00042, 2.772e-02, 0.0f, 0.0f, 0.0f}); // 150 mum Be + layers.push_back(DetLayer{"ddd0", 0.5, 250, x0IT, xrhoIB, resRPhiIT, resZIT, eff}); + layers.push_back(DetLayer{"ddd1", 1.2, 250, x0IT, xrhoIB, resRPhiIT, resZIT, eff}); + layers.push_back(DetLayer{"ddd2", 2.5, 250, x0IT, xrhoIB, resRPhiIT, resZIT, eff}); + layers.push_back(DetLayer{"bpipe1", 5.7, 250, 0.0014, 9.24e-02, 0.0f, 0.0f, 0.0f}); // 500 mum Be + layers.push_back(DetLayer{"ddd3", 7., 250, x0OT, xrhoOB, resRPhiOT, resZOT, eff}); + layers.push_back(DetLayer{"ddd4", 10., 250, x0OT, xrhoOB, resRPhiOT, resZOT, eff}); + layers.push_back(DetLayer{"ddd5", 13., 250, x0OT, xrhoOB, resRPhiOT, resZOT, eff}); + layers.push_back(DetLayer{"ddd6", 16., 250, x0OT, xrhoOB, resRPhiOT, resZOT, eff}); + layers.push_back(DetLayer{"ddd7", 25., 250, x0OT, xrhoOB, resRPhiOT, resZOT, eff}); + layers.push_back(DetLayer{"ddd8", 40., 250, x0OT, xrhoOB, resRPhiOT, resZOT, eff}); + layers.push_back(DetLayer{"ddd9", 45., 250, x0OT, xrhoOB, resRPhiOT, resZOT, eff}); +} + +void FastTracker::AddSiliconALICE3v1() +{ + LOG(info) << " Adding ALICE 3 v1 ITS layers"; + float x0IT = 0.001; // 0.1% + float x0OT = 0.005; // 0.5% + float xrhoIB = 2.3292e-02; // 100 mum Si + float xrhoOB = 2.3292e-01; // 1000 mum Si + + float resRPhiIT = 0.00025; // 2.5 mum + float resZIT = 0.00025; // 2.5 mum + float resRPhiOT = 0.00100; // 5 mum + float resZOT = 0.00100; // 5 mum + float eff = 1.00; + + layers.push_back(DetLayer{"bpipe0", 0.48, 250, 0.00042, 2.772e-02, 0.0f, 0.0f, 0.0f}); // 150 mum Be + layers.push_back(DetLayer{"B00", 0.5, 250, x0IT, xrhoIB, resRPhiIT, resZIT, eff}); + layers.push_back(DetLayer{"B01", 1.2, 250, x0IT, xrhoIB, resRPhiIT, resZIT, eff}); + layers.push_back(DetLayer{"B02", 2.5, 250, x0IT, xrhoIB, resRPhiIT, resZIT, eff}); + layers.push_back(DetLayer{"bpipe1", 3.7, 250, 0.0014, 9.24e-02, 0.0f, 0.0f, 0.0f}); // 500 mum Be + layers.push_back(DetLayer{"B03", 3.75, 250, x0OT, xrhoOB, resRPhiOT, resZOT, eff}); + layers.push_back(DetLayer{"B04", 7., 250, x0OT, xrhoOB, resRPhiOT, resZOT, eff}); + layers.push_back(DetLayer{"B05", 12., 250, x0OT, xrhoOB, resRPhiOT, resZOT, eff}); + layers.push_back(DetLayer{"B06", 20., 250, x0OT, xrhoOB, resRPhiOT, resZOT, eff}); + layers.push_back(DetLayer{"B07", 30., 250, x0OT, xrhoOB, resRPhiOT, resZOT, eff}); + layers.push_back(DetLayer{"B08", 45., 250, x0OT, xrhoOB, resRPhiOT, resZOT, eff}); + layers.push_back(DetLayer{"B09", 60., 250, x0OT, xrhoOB, resRPhiOT, resZOT, eff}); + layers.push_back(DetLayer{"B10", 80., 250, x0OT, xrhoOB, resRPhiOT, resZOT, eff}); +} + +void FastTracker::AddTPC(float phiResMean, float zResMean) +{ + LOG(info) << " Adding standard time projection chamber"; + + // porting of DetectorK::AddTPC + // see here: + // https://github.com/AliceO2Group/DelphesO2/blob/master/src/DetectorK/DetectorK.cxx + // % Radiation Lengths ... Average per TPC row (i.e. total/159 ) + const int kNPassiveBound = 2; + const float radLBoundary[kNPassiveBound] = {1.692612e-01, 8.711904e-02}; + const float xrhoBoundary[kNPassiveBound] = {6.795774e+00, 3.111401e+00}; + const float rBoundary[kNPassiveBound] = {50, 70.0}; // cm + + float radLPerRow = 0.000036; + + float tpcInnerRadialPitch = 0.75; // cm + float tpcMiddleRadialPitch = 1.0; // cm + float tpcOuterRadialPitch = 1.5; // cm + float innerRows = 63; + float middleRows = 64; + float outerRows = 32; + float tpcRows = (innerRows + middleRows + outerRows); + float rowOneRadius = 85.2; // cm + float row64Radius = 135.1; // cm + float row128Radius = 199.2; // cm + + float zLength = 250.0f; // to be checked + + // add boundaries between ITS and TPC + for (int i = 0; i < kNPassiveBound; i++) { + AddLayer(Form("tpc_boundary%d", i), rBoundary[i], zLength, radLBoundary[i], xrhoBoundary[i]); // dummy errors + } + for (Int_t k = 0; k < tpcRows; k++) { + Float_t rowRadius = 0; + if (k < innerRows) + rowRadius = rowOneRadius + k * tpcInnerRadialPitch; + else if (k >= innerRows && k < (innerRows + middleRows)) + rowRadius = row64Radius + (k - innerRows + 1) * tpcMiddleRadialPitch; + else if (k >= (innerRows + middleRows) && k < tpcRows) + rowRadius = row128Radius + (k - innerRows - middleRows + 1) * tpcOuterRadialPitch; + + AddLayer(Form("tpc_%d", k), rowRadius, zLength, radLPerRow, 0, phiResMean, zResMean, 1.0f); + } +} + +// function to provide a reconstructed track from a perfect input track +// returns number of intercepts (generic for now) +int FastTracker::FastTrack(o2::track::TrackParCov inputTrack, o2::track::TrackParCov& outputTrack) +{ + hits.clear(); + int nIntercepts = 0; + std::array posIni; // provision for != PV + inputTrack.getXYZGlo(posIni); + float initialRadius = std::hypot(posIni[0], posIni[1]); + + // +-~-<*>-~-+-~-<*>-~-+-~-<*>-~-+-~-<*>-~-+-~-<*>-~-+ + // Outward pass to find intercepts + int firstLayerReached = -1; + int lastLayerReached = -1; + for (uint32_t il = 0; il < layers.size(); il++) { + // check if layer is doable + if (layers[il].r < initialRadius) + continue; // this layer should not be attempted, but go ahead + if (layers[il].eff < 1e-5) + continue; // inert layer, skip + + // check if layer is reached + float targetX = 1e+3; + inputTrack.getXatLabR(layers[il].r, targetX, magneticField); + if (targetX > 999) + break; // failed to find intercept + + if (!inputTrack.propagateTo(targetX, magneticField)) { + break; // failed to propagate + } + if (std::abs(inputTrack.getZ()) > layers[il].z && applyZacceptance) { + break; // out of acceptance bounds + } + + // layer is reached + if (firstLayerReached < 0) + firstLayerReached = il; + lastLayerReached = il; + nIntercepts++; + } + + // +-~-<*>-~-+-~-<*>-~-+-~-<*>-~-+-~-<*>-~-+-~-<*>-~-+ + // initialize track at outer point + new (&outputTrack)(o2::track::TrackParCov)(inputTrack); + + // Enlarge covariance matrix + std::array trPars = {0.}; + for (int ip = 0; ip < 5; ip++) { + trPars[ip] = outputTrack.getParam(ip); + } + std::array largeCov = {0.}; + enum { kY, + kZ, + kSnp, + kTgl, + kPtI }; // track parameter aliases + enum { kY2, + kYZ, + kZ2, + kYSnp, + kZSnp, + kSnp2, + kYTgl, + kZTgl, + kSnpTgl, + kTgl2, + kYPtI, + kZPtI, + kSnpPtI, + kTglPtI, + kPtI2 }; // cov.matrix aliases + const double kLargeErr2Coord = 5 * 5; + const double kLargeErr2Dir = 0.7 * 0.7; + const double kLargeErr2PtI = 30.5 * 30.5; + for (int ic = 15; ic--;) + largeCov[ic] = 0.; + largeCov[kY2] = largeCov[kZ2] = kLargeErr2Coord; + largeCov[kSnp2] = largeCov[kTgl2] = kLargeErr2Dir; + largeCov[kPtI2] = kLargeErr2PtI * trPars[kPtI] * trPars[kPtI]; + + outputTrack.setCov(largeCov); + outputTrack.checkCovariance(); + + // +-~-<*>-~-+-~-<*>-~-+-~-<*>-~-+-~-<*>-~-+-~-<*>-~-+ + // Inward pass to calculate covariances + for (int il = lastLayerReached; il >= firstLayerReached; il--) { + float targetX = 1e+3; + inputTrack.getXatLabR(layers[il].r, targetX, magneticField); + if (targetX > 999) + continue; // failed to find intercept + + if (!inputTrack.propagateTo(targetX, magneticField)) { + continue; // failed to propagate + } + if (std::abs(inputTrack.getZ()) > layers[il].z && applyZacceptance) { + continue; // out of acceptance bounds but continue inwards + } + + // get perfect data point position + std::array spacePoint; + inputTrack.getXYZGlo(spacePoint); + std::vector thisHit = {spacePoint[0], spacePoint[1], spacePoint[2]}; + + // towards adding cluster: move to track alpha + double alpha = outputTrack.getAlpha(); + double xyz1[3]{ + TMath::Cos(alpha) * spacePoint[0] + TMath::Sin(alpha) * spacePoint[1], + -TMath::Sin(alpha) * spacePoint[0] + TMath::Cos(alpha) * spacePoint[1], + spacePoint[2]}; + if (!(outputTrack.propagateTo(xyz1[0], magneticField))) + continue; + + const o2::track::TrackParametrization::dim2_t hitpoint = { + static_cast(xyz1[1]), + static_cast(xyz1[2])}; + const o2::track::TrackParametrization::dim3_t hitpointcov = {layers[il].resRPhi * layers[il].resRPhi, 0.f, layers[il].resZ * layers[il].resZ}; + outputTrack.update(hitpoint, hitpointcov); + + hits.push_back(thisHit); + outputTrack.checkCovariance(); + } + + // backpropagate to original radius + float finalX = 1e+3; + outputTrack.getXatLabR(initialRadius, finalX, magneticField); + if (finalX > 999) + return -3; // failed to find intercept + + if (!outputTrack.propagateTo(finalX, magneticField)) { + return -4; // failed to propagate + } + + // only attempt to continue if intercepts are at least four + if (nIntercepts < 4) + return nIntercepts; + + // Use covariance matrix based smearing + std::array covMat = {0.}; + for (int ii = 0; ii < 15; ii++) + covMat[ii] = outputTrack.getCov()[ii]; + TMatrixDSym m(5); + double fcovm[5][5]; + + for (int ii = 0, k = 0; ii < 5; ++ii) { + for (int j = 0; j < ii + 1; ++j, ++k) { + fcovm[ii][j] = covMat[k]; + fcovm[j][ii] = covMat[k]; + } + } + + // evaluate ruben's conditional, regularise + bool makePositiveDefinite = (covMatFactor > -1e-5); // apply fix + bool rubenConditional = false; + for (int ii = 0; ii < 5; ii++) { + for (int jj = 0; jj < 5; jj++) { + if (ii == jj) + continue; // don't evaluate diagonals + if (fcovm[ii][jj] * fcovm[ii][jj] > std::abs(fcovm[ii][ii] * fcovm[jj][jj])) { + rubenConditional = true; + if (makePositiveDefinite) { + fcovm[ii][jj] = TMath::Sign(1, fcovm[ii][jj]) * covMatFactor * sqrt(std::abs(fcovm[ii][ii] * fcovm[jj][jj])); + } + } + } + } + + // Should have a valid cov matrix now + m.SetMatrixArray(reinterpret_cast(fcovm)); + TMatrixDSymEigen eigen(m); + TMatrixD eigVec = eigen.GetEigenVectors(); + TVectorD eigVal = eigen.GetEigenValues(); + bool negEigVal = false; + for (int ii = 0; ii < 5; ii++) { + if (eigVal[ii] < 0.0f) + negEigVal = true; + } + + if (negEigVal && rubenConditional && makePositiveDefinite) { + if (verboseLevel > 0) { + LOG(info) << "WARNING: this diagonalization (at pt = " << inputTrack.getPt() << ") has negative eigenvalues despite Ruben's fix! Please be careful!"; + LOG(info) << "Printing info:"; + LOG(info) << "Kalman updates: " << nIntercepts; + LOG(info) << "Cov matrix: "; + m.Print(); + } + covMatNotOK++; + nIntercepts = -1; // mark as problematic so that it isn't used + return -1; + } + covMatOK++; + + // transform parameter vector and smear + double params_[5]; + for (int ii = 0; ii < 5; ++ii) { + double val = 0.; + for (int j = 0; j < 5; ++j) + val += eigVec[j][ii] * outputTrack.getParam(j); + // smear parameters according to eigenvalues + params_[ii] = gRandom->Gaus(val, sqrt(eigVal[ii])); + } + + // invert eigenvector matrix + eigVec.Invert(); + // transform back params vector + for (int ii = 0; ii < 5; ++ii) { + double val = 0.; + for (int j = 0; j < 5; ++j) + val += eigVec[j][ii] * params_[j]; + outputTrack.setParam(val, ii); + } + // should make a sanity check that par[2] sin(phi) is in [-1, 1] + if (fabs(outputTrack.getParam(2)) > 1.) { + LOG(info) << " --- smearTrack failed sin(phi) sanity check: " << outputTrack.getParam(2); + return -2; + } + + return nIntercepts; +} +// +-~-<*>-~-+-~-<*>-~-+-~-<*>-~-+-~-<*>-~-+-~-<*>-~-+-~-<*>-~-+-~-<*>-~-+-~-<*>-~-+ + +} /* namespace fastsim */ +} /* namespace o2 */ + +ClassImp(o2::fastsim::FastTracker); diff --git a/ALICE3/Core/FastTracker.h b/ALICE3/Core/FastTracker.h new file mode 100644 index 00000000000..e8ee7f500e6 --- /dev/null +++ b/ALICE3/Core/FastTracker.h @@ -0,0 +1,67 @@ +// Copyright 2019-2020 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. + +#ifndef ALICE3_CORE_FASTTRACKER_H_ +#define ALICE3_CORE_FASTTRACKER_H_ + +#include // not a system header but megalinter thinks so +#include +#include "DetLayer.h" +#include "ReconstructionDataFormats/Track.h" + +namespace o2 +{ +namespace fastsim +{ + +// +-~-<*>-~-+-~-<*>-~-+-~-<*>-~-+-~-<*>-~-+-~-<*>-~-+-~-<*>-~-+-~-<*>-~-+-~-<*>-~-+ + +// this class implements a synthetic smearer that allows +// for on-demand smearing of TrackParCovs in a certain flexible t +// detector layout. +class FastTracker +{ + public: + // Constructor/destructor + FastTracker(); + virtual ~FastTracker() {} + + void AddLayer(TString name, float r, float z, float x0, float xrho, float resRPhi = 0.0f, float resZ = 0.0f, float eff = 0.0f); + + void AddSiliconALICE3v4(); + void AddSiliconALICE3v1(); + void AddTPC(float phiResMean, float zResMean); + + void Print(); + int FastTrack(o2::track::TrackParCov inputTrack, o2::track::TrackParCov& outputTrack); + + // Definition of detector layers + std::vector layers; + std::vector> hits; // bookkeep last added hits + + // operational + float magneticField; // in kiloGauss (5 = 0.5T, etc) + bool applyZacceptance; // check z acceptance or not + float covMatFactor; // covmat off-diagonal factor to use for covmat fix (negative: no factor) + int verboseLevel; // 0: not verbose, >0 more verbose + + uint64_t covMatOK; // cov mat has negative eigenvals + uint64_t covMatNotOK; // cov mat has negative eigenvals + + ClassDef(FastTracker, 1); +}; + +// +-~-<*>-~-+-~-<*>-~-+-~-<*>-~-+-~-<*>-~-+-~-<*>-~-+-~-<*>-~-+-~-<*>-~-+-~-<*>-~-+ + +} // namespace fastsim +} // namespace o2 + +#endif // ALICE3_CORE_FASTTRACKER_H_ diff --git a/ALICE3/Core/FastTrackerLinkDef.h b/ALICE3/Core/FastTrackerLinkDef.h new file mode 100644 index 00000000000..a69755b7e92 --- /dev/null +++ b/ALICE3/Core/FastTrackerLinkDef.h @@ -0,0 +1,21 @@ +// Copyright 2019-2020 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. + +#ifndef ALICE3_CORE_FASTTRACKERLINKDEF_H_ +#define ALICE3_CORE_FASTTRACKERLINKDEF_H_ + +#pragma link off all globals; +#pragma link off all classes; +#pragma link off all functions; + +#pragma link C++ class o2::fastsim::FastTracker + ; + +#endif // ALICE3_CORE_FASTTRACKERLINKDEF_H_ diff --git a/ALICE3/TableProducer/OTF/CMakeLists.txt b/ALICE3/TableProducer/OTF/CMakeLists.txt index 55427691f4f..d20b6c6ba5c 100644 --- a/ALICE3/TableProducer/OTF/CMakeLists.txt +++ b/ALICE3/TableProducer/OTF/CMakeLists.txt @@ -11,7 +11,7 @@ o2physics_add_dpl_workflow(onthefly-tracker SOURCES onTheFlyTracker.cxx - PUBLIC_LINK_LIBRARIES O2::Framework O2::DetectorsBase O2Physics::AnalysisCore O2::ReconstructionDataFormats O2::DetectorsCommonDataFormats O2::DetectorsVertexing O2::DCAFitter O2Physics::ALICE3Core + PUBLIC_LINK_LIBRARIES O2::Framework O2::DetectorsBase O2Physics::AnalysisCore O2::ReconstructionDataFormats O2::DetectorsCommonDataFormats O2::DetectorsVertexing O2::DCAFitter O2Physics::ALICE3Core O2Physics::FastTracker COMPONENT_NAME Analysis) o2physics_add_dpl_workflow(onthefly-tofpid diff --git a/ALICE3/TableProducer/OTF/onTheFlyTracker.cxx b/ALICE3/TableProducer/OTF/onTheFlyTracker.cxx index 11f2746fcf1..9ee71307610 100644 --- a/ALICE3/TableProducer/OTF/onTheFlyTracker.cxx +++ b/ALICE3/TableProducer/OTF/onTheFlyTracker.cxx @@ -56,6 +56,7 @@ #include "ITStracking/VertexerTraits.h" #include "ALICE3/Core/DelphesO2TrackSmearer.h" +#include "ALICE3/Core/FastTracker.h" #include "ALICE3/DataModel/collisionAlice3.h" #include "ALICE3/DataModel/tracksAlice3.h" #include "ALICE3/DataModel/OTFStrangeness.h" @@ -118,37 +119,40 @@ struct OnTheFlyTracker { Configurable lutTr{"lutTr", "lutCovm.tr.dat", "LUT for tritons"}; Configurable lutHe3{"lutHe3", "lutCovm.he3.dat", "LUT for Helium-3"}; - Configurable lutPi0{"lutPi0", "lutCovm.pi.20kG.rmin20.geometry_v0.dat", "LUT for pions without layer 0"}; - Configurable lutPi1{"lutPi1", "lutCovm.pi.20kG.rmin20.geometry_v1.dat", "LUT for pions without layer 1"}; - Configurable lutPi2{"lutPi2", "lutCovm.pi.20kG.rmin20.geometry_v2.dat", "LUT for pions without layer 2"}; - Configurable lutPi3{"lutPi3", "lutCovm.pi.20kG.rmin20.geometry_v3.dat", "LUT for pions without layer 3"}; - Configurable lutPi4{"lutPi4", "lutCovm.pi.20kG.rmin20.geometry_v4.dat", "LUT for pions without layer 4"}; - Configurable lutPi5{"lutPi5", "lutCovm.pi.20kG.rmin20.geometry_v5.dat", "LUT for pions without layer 5"}; - Configurable lutPr0{"lutPr0", "lutCovm.pr.20kG.rmin20.geometry_v0.dat", "LUT for protons without layer 0"}; - Configurable lutPr1{"lutPr1", "lutCovm.pr.20kG.rmin20.geometry_v1.dat", "LUT for protons without layer 1"}; - Configurable lutPr2{"lutPr2", "lutCovm.pr.20kG.rmin20.geometry_v2.dat", "LUT for protons without layer 2"}; - Configurable lutPr3{"lutPr3", "lutCovm.pr.20kG.rmin20.geometry_v3.dat", "LUT for protons without layer 3"}; - Configurable lutPr4{"lutPr4", "lutCovm.pr.20kG.rmin20.geometry_v4.dat", "LUT for protons without layer 4"}; - Configurable lutPr5{"lutPr5", "lutCovm.pr.20kG.rmin20.geometry_v5.dat", "LUT for protons without layer 5"}; - - ConfigurableAxis axisMomentum{"axisMomentum", {VARIABLE_WIDTH, 0.0f, 0.1f, 0.2f, 0.3f, 0.4f, 0.5f, 0.6f, 0.7f, 0.8f, 0.9f, 1.0f, 1.1f, 1.2f, 1.3f, 1.4f, 1.5f, 1.6f, 1.7f, 1.8f, 1.9f, 2.0f, 2.2f, 2.4f, 2.6f, 2.8f, 3.0f, 3.2f, 3.4f, 3.6f, 3.8f, 4.0f, 4.4f, 4.8f, 5.2f, 5.6f, 6.0f, 6.5f, 7.0f, 7.5f, 8.0f, 9.0f, 10.0f, 11.0f, 12.0f, 13.0f, 14.0f, 15.0f, 17.0f, 19.0f, 21.0f, 23.0f, 25.0f, 30.0f, 35.0f, 40.0f, 50.0f}, "#it{p} (GeV/#it{c})"}; - ConfigurableAxis axisNVertices{"axisNVertices", {20, -0.5, 19.5}, "N_{vertices}"}; - ConfigurableAxis axisMultiplicity{"axisMultiplicity", {100, -0.5, 99.5}, "N_{contributors}"}; - ConfigurableAxis axisVertexZ{"axisVertexZ", {40, -20, 20}, "vertex Z (cm)"}; - ConfigurableAxis axisDCA{"axisDCA", {400, -200, 200}, "DCA (#mum)"}; - ConfigurableAxis axisX{"axisX", {250, -50, 200}, "track X (cm)"}; - ConfigurableAxis axisRadius{"axisRadius", {55, 0.01, 100}, "decay radius"}; - ConfigurableAxis axisLambdaMass{"axisLambdaMass", {200, 1.101f, 1.131f}, ""}; - ConfigurableAxis axisXiMass{"axisXiMass", {200, 1.22f, 1.42f}, ""}; - - ConfigurableAxis axisDeltaPt{"axisDeltaPt", {200, -1.0f, +1.0f}, "#Delta p_{T}"}; - ConfigurableAxis axisDeltaEta{"axisDeltaEta", {200, -0.5f, +0.5f}, "#Delta #eta"}; + struct : ConfigurableGroup { + ConfigurableAxis axisMomentum{"axisMomentum", {VARIABLE_WIDTH, 0.0f, 0.1f, 0.2f, 0.3f, 0.4f, 0.5f, 0.6f, 0.7f, 0.8f, 0.9f, 1.0f, 1.1f, 1.2f, 1.3f, 1.4f, 1.5f, 1.6f, 1.7f, 1.8f, 1.9f, 2.0f, 2.2f, 2.4f, 2.6f, 2.8f, 3.0f, 3.2f, 3.4f, 3.6f, 3.8f, 4.0f, 4.4f, 4.8f, 5.2f, 5.6f, 6.0f, 6.5f, 7.0f, 7.5f, 8.0f, 9.0f, 10.0f, 11.0f, 12.0f, 13.0f, 14.0f, 15.0f, 17.0f, 19.0f, 21.0f, 23.0f, 25.0f, 30.0f, 35.0f, 40.0f, 50.0f}, "#it{p} (GeV/#it{c})"}; + ConfigurableAxis axisNVertices{"axisNVertices", {20, -0.5, 19.5}, "N_{vertices}"}; + ConfigurableAxis axisMultiplicity{"axisMultiplicity", {100, -0.5, 99.5}, "N_{contributors}"}; + ConfigurableAxis axisVertexZ{"axisVertexZ", {40, -20, 20}, "vertex Z (cm)"}; + ConfigurableAxis axisDCA{"axisDCA", {400, -200, 200}, "DCA (#mum)"}; + ConfigurableAxis axisX{"axisX", {250, -50, 200}, "track X (cm)"}; + ConfigurableAxis axisDecayRadius{"axisDecayRadius", {55, 0.01, 100}, "decay radius"}; + ConfigurableAxis axisLambdaMass{"axisLambdaMass", {200, 1.101f, 1.131f}, ""}; + ConfigurableAxis axisXiMass{"axisXiMass", {200, 1.22f, 1.42f}, ""}; + + ConfigurableAxis axisDeltaPt{"axisDeltaPt", {200, -1.0f, +1.0f}, "#Delta p_{T}"}; + ConfigurableAxis axisDeltaEta{"axisDeltaEta", {200, -0.5f, +0.5f}, "#Delta #eta"}; + + ConfigurableAxis axisRadius{"axisRadius", {2500, 0.0f, +250.0f}, "R (cm)"}; + ConfigurableAxis axisZ{"axisZ", {100, -250.0f, +250.0f}, "Z (cm)"}; + } axes; + + // for topo var QA + struct : ConfigurableGroup { + std::string prefix = "fastTrackerSettings"; // JSON group name + Configurable minSiliconHits{"minSiliconHits", 4, "minimum number of silicon hits to accept track"}; + Configurable alice3detector{"alice3detector", 0, "0: ALICE 3 v1, 1: ALICE 3 v4"}; + Configurable applyZacceptance{"applyZacceptance", false, "apply z limits to detector layers or not"}; + } fastTrackerSettings; // allows for gap between peak and bg in case someone wants to using PVertex = o2::dataformats::PrimaryVertex; // for secondary vertex finding o2::vertexing::DCAFitterN<2> fitter; + // FastTracker machinery + o2::fastsim::FastTracker fastTracker; + // Class to hold the track information for the O2 vertexing class TrackAlice3 : public o2::track::TrackParCov { @@ -264,103 +268,6 @@ struct OnTheFlyTracker { // smear un-reco'ed tracks if asked to do so mSmearer.skipUnreconstructed(static_cast(!processUnreconstructedTracks)); - - if (treatXi) { - std::map mapPdgLut0; - const char* lutPiChar0 = lutPi0->c_str(); - const char* lutPrChar0 = lutPr0->c_str(); - LOGF(info, "Load more pion lut files .....: %s", lutPiChar0); - LOGF(info, "Load more proton lut files ...: %s", lutPrChar0); - mapPdgLut0.insert(std::make_pair(211, lutPiChar0)); - mapPdgLut0.insert(std::make_pair(2212, lutPrChar0)); - - std::map mapPdgLut1; - const char* lutPiChar1 = lutPi1->c_str(); - const char* lutPrChar1 = lutPr1->c_str(); - LOGF(info, "Load more pion lut files .....: %s", lutPiChar1); - LOGF(info, "Load more proton lut files ...: %s", lutPrChar1); - mapPdgLut1.insert(std::make_pair(211, lutPiChar1)); - mapPdgLut1.insert(std::make_pair(2212, lutPrChar1)); - - std::map mapPdgLut2; - const char* lutPiChar2 = lutPi2->c_str(); - const char* lutPrChar2 = lutPr2->c_str(); - LOGF(info, "Load more pion lut files .....: %s", lutPiChar2); - LOGF(info, "Load more proton lut files ...: %s", lutPrChar2); - mapPdgLut2.insert(std::make_pair(211, lutPiChar2)); - mapPdgLut2.insert(std::make_pair(2212, lutPrChar2)); - - std::map mapPdgLut3; - const char* lutPiChar3 = lutPi3->c_str(); - const char* lutPrChar3 = lutPr3->c_str(); - LOGF(info, "Load more pion lut files .....: %s", lutPiChar3); - LOGF(info, "Load more proton lut files ...: %s", lutPrChar3); - mapPdgLut3.insert(std::make_pair(211, lutPiChar3)); - mapPdgLut3.insert(std::make_pair(2212, lutPrChar3)); - - std::map mapPdgLut4; - const char* lutPiChar4 = lutPi4->c_str(); - const char* lutPrChar4 = lutPr4->c_str(); - LOGF(info, "Load more pion lut files .....: %s", lutPiChar4); - LOGF(info, "Load more proton lut files ...: %s", lutPrChar4); - mapPdgLut4.insert(std::make_pair(211, lutPiChar4)); - mapPdgLut4.insert(std::make_pair(2212, lutPrChar4)); - - std::map mapPdgLut5; - const char* lutPiChar5 = lutPi5->c_str(); - const char* lutPrChar5 = lutPr5->c_str(); - LOGF(info, "Load more pion lut files .....: %s", lutPiChar5); - LOGF(info, "Load more proton lut files ...: %s", lutPrChar5); - mapPdgLut5.insert(std::make_pair(211, lutPiChar5)); - mapPdgLut5.insert(std::make_pair(2212, lutPrChar5)); - - for (auto e : mapPdgLut0) { - if (!mSmearer0.loadTable(e.first, e.second)) { - LOG(fatal) << "Having issue with loading the LUT " << e.first << " " << e.second; - } - } - for (auto e : mapPdgLut1) { - if (!mSmearer1.loadTable(e.first, e.second)) { - LOG(fatal) << "Having issue with loading the LUT " << e.first << " " << e.second; - } - } - for (auto e : mapPdgLut2) { - if (!mSmearer2.loadTable(e.first, e.second)) { - LOG(fatal) << "Having issue with loading the LUT " << e.first << " " << e.second; - } - } - for (auto e : mapPdgLut3) { - if (!mSmearer3.loadTable(e.first, e.second)) { - LOG(fatal) << "Having issue with loading the LUT " << e.first << " " << e.second; - } - } - for (auto e : mapPdgLut4) { - if (!mSmearer4.loadTable(e.first, e.second)) { - LOG(fatal) << "Having issue with loading the LUT " << e.first << " " << e.second; - } - } - for (auto e : mapPdgLut5) { - if (!mSmearer5.loadTable(e.first, e.second)) { - LOG(fatal) << "Having issue with loading the LUT " << e.first << " " << e.second; - } - } - - // interpolate efficiencies if requested to do so - mSmearer0.interpolateEfficiency(static_cast(interpolateLutEfficiencyVsNch)); - mSmearer1.interpolateEfficiency(static_cast(interpolateLutEfficiencyVsNch)); - mSmearer2.interpolateEfficiency(static_cast(interpolateLutEfficiencyVsNch)); - mSmearer3.interpolateEfficiency(static_cast(interpolateLutEfficiencyVsNch)); - mSmearer4.interpolateEfficiency(static_cast(interpolateLutEfficiencyVsNch)); - mSmearer5.interpolateEfficiency(static_cast(interpolateLutEfficiencyVsNch)); - - // smear un-reco'ed tracks if asked to do so - mSmearer0.skipUnreconstructed(static_cast(!processUnreconstructedTracks)); - mSmearer1.skipUnreconstructed(static_cast(!processUnreconstructedTracks)); - mSmearer2.skipUnreconstructed(static_cast(!processUnreconstructedTracks)); - mSmearer3.skipUnreconstructed(static_cast(!processUnreconstructedTracks)); - mSmearer4.skipUnreconstructed(static_cast(!processUnreconstructedTracks)); - mSmearer5.skipUnreconstructed(static_cast(!processUnreconstructedTracks)); - } } // Basic QA @@ -374,60 +281,66 @@ struct OnTheFlyTracker { hNaN->GetYaxis()->SetBinLabel(1, "Smear NaN"); hNaN->GetYaxis()->SetBinLabel(2, "Smear OK"); - histos.add("hPtGenerated", "hPtGenerated", kTH1F, {axisMomentum}); - histos.add("hPtGeneratedEl", "hPtGeneratedEl", kTH1F, {axisMomentum}); - histos.add("hPtGeneratedPi", "hPtGeneratedPi", kTH1F, {axisMomentum}); - histos.add("hPtGeneratedKa", "hPtGeneratedKa", kTH1F, {axisMomentum}); - histos.add("hPtGeneratedPr", "hPtGeneratedPr", kTH1F, {axisMomentum}); - histos.add("hPtReconstructed", "hPtReconstructed", kTH1F, {axisMomentum}); - histos.add("hPtReconstructedEl", "hPtReconstructedEl", kTH1F, {axisMomentum}); - histos.add("hPtReconstructedPi", "hPtReconstructedPi", kTH1F, {axisMomentum}); - histos.add("hPtReconstructedKa", "hPtReconstructedKa", kTH1F, {axisMomentum}); - histos.add("hPtReconstructedPr", "hPtReconstructedPr", kTH1F, {axisMomentum}); + auto hCovMatOK = histos.add("hCovMatOK", "hCovMatOK", kTH1D, {{2, -0.5f, 1.5f}}); + hCovMatOK->GetXaxis()->SetBinLabel(1, "Not OK"); + hCovMatOK->GetXaxis()->SetBinLabel(2, "OK"); + + histos.add("hPtGenerated", "hPtGenerated", kTH1F, {axes.axisMomentum}); + histos.add("hPtGeneratedEl", "hPtGeneratedEl", kTH1F, {axes.axisMomentum}); + histos.add("hPtGeneratedPi", "hPtGeneratedPi", kTH1F, {axes.axisMomentum}); + histos.add("hPtGeneratedKa", "hPtGeneratedKa", kTH1F, {axes.axisMomentum}); + histos.add("hPtGeneratedPr", "hPtGeneratedPr", kTH1F, {axes.axisMomentum}); + histos.add("hPtReconstructed", "hPtReconstructed", kTH1F, {axes.axisMomentum}); + histos.add("hPtReconstructedEl", "hPtReconstructedEl", kTH1F, {axes.axisMomentum}); + histos.add("hPtReconstructedPi", "hPtReconstructedPi", kTH1F, {axes.axisMomentum}); + histos.add("hPtReconstructedKa", "hPtReconstructedKa", kTH1F, {axes.axisMomentum}); + histos.add("hPtReconstructedPr", "hPtReconstructedPr", kTH1F, {axes.axisMomentum}); // Collision QA - histos.add("hPVz", "hPVz", kTH1F, {axisVertexZ}); - histos.add("hLUTMultiplicity", "hLUTMultiplicity", kTH1F, {axisMultiplicity}); - histos.add("hSimMultiplicity", "hSimMultiplicity", kTH1F, {axisMultiplicity}); - histos.add("hRecoMultiplicity", "hRecoMultiplicity", kTH1F, {axisMultiplicity}); + histos.add("hPVz", "hPVz", kTH1F, {axes.axisVertexZ}); + histos.add("hLUTMultiplicity", "hLUTMultiplicity", kTH1F, {axes.axisMultiplicity}); + histos.add("hSimMultiplicity", "hSimMultiplicity", kTH1F, {axes.axisMultiplicity}); + histos.add("hRecoMultiplicity", "hRecoMultiplicity", kTH1F, {axes.axisMultiplicity}); if (doExtraQA) { - histos.add("h2dVerticesVsContributors", "h2dVerticesVsContributors", kTH2F, {axisMultiplicity, axisNVertices}); - histos.add("hRecoVsSimMultiplicity", "hRecoVsSimMultiplicity", kTH2F, {axisMultiplicity, axisMultiplicity}); - histos.add("h2dDCAxy", "h2dDCAxy", kTH2F, {axisMomentum, axisDCA}); + histos.add("h2dVerticesVsContributors", "h2dVerticesVsContributors", kTH2F, {axes.axisMultiplicity, axes.axisNVertices}); + histos.add("hRecoVsSimMultiplicity", "hRecoVsSimMultiplicity", kTH2F, {axes.axisMultiplicity, axes.axisMultiplicity}); + histos.add("h2dDCAxy", "h2dDCAxy", kTH2F, {axes.axisMomentum, axes.axisDCA}); - histos.add("hSimTrackX", "hSimTrackX", kTH1F, {axisX}); - histos.add("hRecoTrackX", "hRecoTrackX", kTH1F, {axisX}); - histos.add("hTrackXatDCA", "hTrackXatDCA", kTH1F, {axisX}); + histos.add("hSimTrackX", "hSimTrackX", kTH1F, {axes.axisX}); + histos.add("hRecoTrackX", "hRecoTrackX", kTH1F, {axes.axisX}); + histos.add("hTrackXatDCA", "hTrackXatDCA", kTH1F, {axes.axisX}); } if (doXiQA) { histos.add("hXiBuilding", "hXiBuilding", kTH1F, {{10, -0.5f, 9.5f}}); - histos.add("hGenXi", "hGenXi", kTH2F, {axisRadius, axisMomentum}); - histos.add("hRecoXi", "hRecoXi", kTH2F, {axisRadius, axisMomentum}); + histos.add("hGenXi", "hGenXi", kTH2F, {axes.axisDecayRadius, axes.axisMomentum}); + histos.add("hRecoXi", "hRecoXi", kTH2F, {axes.axisDecayRadius, axes.axisMomentum}); - histos.add("hGenPiFromXi", "hGenPiFromXi", kTH2F, {axisRadius, axisMomentum}); - histos.add("hGenPiFromL0", "hGenPiFromL0", kTH2F, {axisRadius, axisMomentum}); - histos.add("hGenPrFromL0", "hGenPrFromL0", kTH2F, {axisRadius, axisMomentum}); - histos.add("hRecoPiFromXi", "hRecoPiFromXi", kTH2F, {axisRadius, axisMomentum}); - histos.add("hRecoPiFromL0", "hRecoPiFromL0", kTH2F, {axisRadius, axisMomentum}); - histos.add("hRecoPrFromL0", "hRecoPrFromL0", kTH2F, {axisRadius, axisMomentum}); + histos.add("hGenPiFromXi", "hGenPiFromXi", kTH2F, {axes.axisDecayRadius, axes.axisMomentum}); + histos.add("hGenPiFromL0", "hGenPiFromL0", kTH2F, {axes.axisDecayRadius, axes.axisMomentum}); + histos.add("hGenPrFromL0", "hGenPrFromL0", kTH2F, {axes.axisDecayRadius, axes.axisMomentum}); + histos.add("hRecoPiFromXi", "hRecoPiFromXi", kTH2F, {axes.axisDecayRadius, axes.axisMomentum}); + histos.add("hRecoPiFromL0", "hRecoPiFromL0", kTH2F, {axes.axisDecayRadius, axes.axisMomentum}); + histos.add("hRecoPrFromL0", "hRecoPrFromL0", kTH2F, {axes.axisDecayRadius, axes.axisMomentum}); // basic mass histograms to see if we're in business - histos.add("hMassLambda", "hMassLambda", kTH1F, {axisLambdaMass}); - histos.add("hMassXi", "hMassXi", kTH1F, {axisXiMass}); + histos.add("hMassLambda", "hMassLambda", kTH1F, {axes.axisLambdaMass}); + histos.add("hMassXi", "hMassXi", kTH1F, {axes.axisXiMass}); // OTF strangeness tracking QA histos.add("hFoundVsFindable", "hFoundVsFindable", kTH2F, {{10, -0.5f, 9.5f}, {10, -0.5f, 9.5f}}); - histos.add("h2dDCAxyCascade", "h2dDCAxyCascade", kTH2F, {axisMomentum, axisDCA}); - histos.add("h2dDCAxyCascadeBachelor", "h2dDCAxyCascadeBachelor", kTH2F, {axisMomentum, axisDCA}); - histos.add("h2dDCAxyCascadeNegative", "h2dDCAxyCascadeNegative", kTH2F, {axisMomentum, axisDCA}); - histos.add("h2dDCAxyCascadePositive", "h2dDCAxyCascadePositive", kTH2F, {axisMomentum, axisDCA}); + histos.add("h2dDCAxyCascade", "h2dDCAxyCascade", kTH2F, {axes.axisMomentum, axes.axisDCA}); + histos.add("h2dDCAxyCascadeBachelor", "h2dDCAxyCascadeBachelor", kTH2F, {axes.axisMomentum, axes.axisDCA}); + histos.add("h2dDCAxyCascadeNegative", "h2dDCAxyCascadeNegative", kTH2F, {axes.axisMomentum, axes.axisDCA}); + histos.add("h2dDCAxyCascadePositive", "h2dDCAxyCascadePositive", kTH2F, {axes.axisMomentum, axes.axisDCA}); + + histos.add("h2dDeltaPtVsPt", "h2dDeltaPtVsPt", kTH2F, {axes.axisMomentum, axes.axisDeltaPt}); + histos.add("h2dDeltaEtaVsPt", "h2dDeltaEtaVsPt", kTH2F, {axes.axisMomentum, axes.axisDeltaEta}); - histos.add("h2dDeltaPtVsPt", "h2dDeltaPtVsPt", kTH2F, {axisMomentum, axisDeltaPt}); - histos.add("h2dDeltaEtaVsPt", "h2dDeltaEtaVsPt", kTH2F, {axisMomentum, axisDeltaEta}); + histos.add("hFastTrackerHits", "hFastTrackerHits", kTH2F, {axes.axisZ, axes.axisRadius}); } LOGF(info, "Initializing magnetic field to value: %.3f kG", static_cast(magneticField)); @@ -475,6 +388,21 @@ struct OnTheFlyTracker { // Set seed for TGenPhaseSpace rand.SetSeed(seed); + + // configure FastTracker + fastTracker.magneticField = magneticField; + fastTracker.applyZacceptance = fastTrackerSettings.applyZacceptance; + + if (fastTrackerSettings.alice3detector == 0) { + fastTracker.AddSiliconALICE3v1(); + } + if (fastTrackerSettings.alice3detector == 1) { + fastTracker.AddSiliconALICE3v4(); + fastTracker.AddTPC(0.1, 0.1); + } + + // print fastTracker settings + fastTracker.Print(); } /// Function to decay the xi @@ -691,8 +619,10 @@ struct OnTheFlyTracker { multiplicityCounter++; const float t = (ir.timeInBCNS + gRandom->Gaus(0., 100.)) * 1e-3; - std::vector xiDaughterTrackParCovs(3); + std::vector xiDaughterTrackParCovsPerfect(3); + std::vector xiDaughterTrackParCovsTracked(3); std::vector isReco(3); + std::vector nHits(3); std::vector smearer = {mSmearer0, mSmearer1, mSmearer2, mSmearer3, mSmearer4, mSmearer5}; if (treatXi && mcParticle.pdgCode() == 3312) { histos.fill(HIST("hXiBuilding"), 0.0f); @@ -700,65 +630,38 @@ struct OnTheFlyTracker { continue; } - histos.fill(HIST("hXiBuilding"), 1.0f); - convertTLorentzVectorToO2Track(-211, decayProducts[0], xiDecayVertex, xiDaughterTrackParCovs[0]); - convertTLorentzVectorToO2Track(-211, decayProducts[1], l0DecayVertex, xiDaughterTrackParCovs[1]); - convertTLorentzVectorToO2Track(2212, decayProducts[2], l0DecayVertex, xiDaughterTrackParCovs[2]); - - // Map daughter to smearer - if (enableSecondarySmearing) { - int firstSmearerIndex = -1; - int secondSmearerIndex = -1; - for (unsigned i = 0; i < layers.size(); i++) { - if (xiDecayRadius2D > layers[i]) { - firstSmearerIndex = i; + convertTLorentzVectorToO2Track(-211, decayProducts[0], xiDecayVertex, xiDaughterTrackParCovsPerfect[0]); + convertTLorentzVectorToO2Track(-211, decayProducts[1], l0DecayVertex, xiDaughterTrackParCovsPerfect[1]); + convertTLorentzVectorToO2Track(2212, decayProducts[2], l0DecayVertex, xiDaughterTrackParCovsPerfect[2]); + + for (int i = 0; i < 3; i++) { + isReco[i] = false; + if (enableSecondarySmearing) { + + nHits[i] = fastTracker.FastTrack(xiDaughterTrackParCovsPerfect[i], xiDaughterTrackParCovsTracked[i]); + + if (nHits[i] >= fastTrackerSettings.minSiliconHits) { + isReco[i] = true; + } else { + continue; // extra sure } - if (l0DecayRadius2D > layers[i]) { - secondSmearerIndex = i; + for (uint32_t ih = 0; ih < fastTracker.hits.size(); ih++) { + histos.fill(HIST("hFastTrackerHits"), fastTracker.hits[ih][2], std::hypot(fastTracker.hits[ih][0], fastTracker.hits[ih][1])); } - } - if (firstSmearerIndex > 5) { - isReco[0] = false; - } else if (firstSmearerIndex == -1) { - isReco[0] = mSmearer.smearTrack(xiDaughterTrackParCovs[0], 211, dNdEta); - } else { - isReco[0] = smearer[firstSmearerIndex].smearTrack(xiDaughterTrackParCovs[0], 211, dNdEta); - } - if (secondSmearerIndex > 5) { - isReco[1] = false; - isReco[2] = false; - } else if (secondSmearerIndex == -1) { - isReco[1] = mSmearer.smearTrack(xiDaughterTrackParCovs[1], 211, dNdEta); - isReco[2] = mSmearer.smearTrack(xiDaughterTrackParCovs[2], 2212, dNdEta); } else { - isReco[1] = smearer[secondSmearerIndex].smearTrack(xiDaughterTrackParCovs[1], 211, dNdEta); - isReco[2] = smearer[secondSmearerIndex].smearTrack(xiDaughterTrackParCovs[2], 2212, dNdEta); + isReco[i] = true; + xiDaughterTrackParCovsTracked[i] = xiDaughterTrackParCovsPerfect[i]; } - } else { - isReco[0] = true; - isReco[1] = true; - isReco[2] = true; - } - for (int i = 0; i < 3; i++) { - if (decayProducts[i].Pt() < minPt) { - isReco[i] = false; - } - if (!isReco[i] && !processUnreconstructedTracks) { - continue; - } - if (TMath::IsNaN(xiDaughterTrackParCovs[i].getZ())) { - histos.fill(HIST("hNaNBookkeeping"), i + 1, 0.0f); - LOGF(info, "Issues with track parametrization %i ! inspect track:", i); - xiDaughterTrackParCovs[i].print(); - isReco[i] = false; // not acceptable + + if (TMath::IsNaN(xiDaughterTrackParCovsTracked[i].getZ())) { continue; } else { histos.fill(HIST("hNaNBookkeeping"), i + 1, 1.0f); } if (isReco[i]) { - tracksAlice3.push_back(TrackAlice3{xiDaughterTrackParCovs[i], mcParticle.globalIndex(), t, 100.f * 1e-3, true, true, i + 2}); + tracksAlice3.push_back(TrackAlice3{xiDaughterTrackParCovsTracked[i], mcParticle.globalIndex(), t, 100.f * 1e-3, true, true, i + 2}); } else { - ghostTracksAlice3.push_back(TrackAlice3{xiDaughterTrackParCovs[i], mcParticle.globalIndex(), t, 100.f * 1e-3, true, true, i + 2}); + ghostTracksAlice3.push_back(TrackAlice3{xiDaughterTrackParCovsTracked[i], mcParticle.globalIndex(), t, 100.f * 1e-3, true, true, i + 2}); } } @@ -793,7 +696,7 @@ struct OnTheFlyTracker { int nCand = 0; bool dcaFitterOK_V0 = true; try { - nCand = fitter.process(xiDaughterTrackParCovs[1], xiDaughterTrackParCovs[2]); + nCand = fitter.process(xiDaughterTrackParCovsTracked[1], xiDaughterTrackParCovsTracked[2]); } catch (...) { // LOG(error) << "Exception caught in DCA fitter process call!"; dcaFitterOK_V0 = false; @@ -845,7 +748,7 @@ struct OnTheFlyTracker { nCand = 0; bool dcaFitterOK_Cascade = true; try { - nCand = fitter.process(v0Track, xiDaughterTrackParCovs[0]); + nCand = fitter.process(v0Track, xiDaughterTrackParCovsTracked[0]); } catch (...) { // LOG(error) << "Exception caught in DCA fitter process call!"; dcaFitterOK_Cascade = false; @@ -1192,7 +1095,12 @@ struct OnTheFlyTracker { cascade.findableClusters, cascade.foundClusters); } - } + + // do bookkeeping of fastTracker tracking + histos.fill(HIST("hCovMatOK"), 0.0f, fastTracker.covMatNotOK); + histos.fill(HIST("hCovMatOK"), 1.0f, fastTracker.covMatOK); + + } // end process }; /// Extends TracksExtra if necessary