diff --git a/ALICE3/Core/FastTracker.cxx b/ALICE3/Core/FastTracker.cxx index 1d06958504b..fe89fcd80fd 100644 --- a/ALICE3/Core/FastTracker.cxx +++ b/ALICE3/Core/FastTracker.cxx @@ -9,13 +9,17 @@ // granted to it by virtue of its status as an Intergovernmental Organization // or submit itself to any jurisdiction. -#include -#include +#include "FastTracker.h" + +#include "ReconstructionDataFormats/TrackParametrization.h" + #include "TMath.h" #include "TMatrixD.h" -#include "TRandom.h" #include "TMatrixDSymEigen.h" -#include "FastTracker.h" +#include "TRandom.h" + +#include +#include namespace o2 { @@ -24,36 +28,6 @@ namespace fastsim // +-~-<*>-~-+-~-<*>-~-+-~-<*>-~-+-~-<*>-~-+-~-<*>-~-+-~-<*>-~-+-~-<*>-~-+-~-<*>-~-+ -FastTracker::FastTracker() -{ - // base constructor - magneticField = 20; // in kiloGauss - applyZacceptance = false; - applyMSCorrection = true; - applyElossCorrection = true; - applyEffCorrection = true; - covMatFactor = 0.99f; - verboseLevel = 0; - - // last fast-tracked track properties - covMatOK = 0; - covMatNotOK = 0; - nIntercepts = 0; - nSiliconPoints = 0; - nGasPoints = 0; - - maxRadiusSlowDet = 10; - integrationTime = 0.02; // ms - crossSectionMinB = 8; - dNdEtaCent = 2200; - dNdEtaMinB = 1; - avgRapidity = 0.45; - sigmaD = 6.0; - luminosity = 1.e27; - otherBackground = 0.0; // [0, 1] - upcBackgroundMultiplier = 1.0; -} - void FastTracker::AddLayer(TString name, float r, float z, float x0, float xrho, float resRPhi, float resZ, float eff, int type) { DetLayer newLayer(name, r, z, x0, xrho, resRPhi, resZ, eff, type); @@ -61,6 +35,17 @@ void FastTracker::AddLayer(TString name, float r, float z, float x0, float xrho, if (newLayer.getEfficiency() > 0.0f && newLayer.isInert()) { LOG(error) << "Layer " << name << " with efficiency > 0.0 should not be inert"; } + // Layers should be ordered by increasing radius, check this + if (!layers.empty() && newLayer.getRadius() < layers.back().getRadius()) { + LOG(fatal) << "Layer " << newLayer << " is not ordered correctly, it should be after layer " << layers.back(); + } + // Layers should all have different names + for (const auto& layer : layers) { + if (layer.getName() == newLayer.getName()) { + LOG(fatal) << "Layer with name " << newLayer.getName() << " already exists in FastTracker layers"; + } + } + // Add the new layer to the layers vector layers.push_back(newLayer); } @@ -78,7 +63,7 @@ DetLayer FastTracker::GetLayer(int layer, bool ignoreBarrelLayers) const return layers[layerIdx]; } -int FastTracker::GetLayerIndex(std::string name) const +int FastTracker::GetLayerIndex(const std::string& name) const { int i = 0; for (const auto& layer : layers) { @@ -211,9 +196,9 @@ float FastTracker::Dist(float z, float r) // https://github.com/AliceO2Group/DelphesO2/blob/master/src/DetectorK/DetectorK.cxx#L743 int index = 1; int nSteps = 301; - double dist = 0.0; - double dz0 = (4 * sigmaD - (-4) * sigmaD / (nSteps = 1)); - double z0 = 0.0; + float dist = 0.0; + float dz0 = (4 * sigmaD - (-4) * sigmaD / (nSteps = 1)); + float z0 = 0.0; for (int i = 0; i < nSteps; i++) { if (i == nSteps - 1) index = 1; @@ -243,7 +228,7 @@ float FastTracker::IntegratedHitDensity(float multiplicity, float radius) // porting of DetektorK::IntegratedHitDensity // see here: // https://github.com/AliceO2Group/DelphesO2/blob/master/src/DetectorK/DetectorK.cxx#L712 - float zdcHz = luminosity * 1.e24 * crossSectionMinB; + float zdcHz = luminosity * 1.e24 * mCrossSectionMinB; float den = zdcHz * integrationTime / 1000. * multiplicity * Dist(0., radius) / (o2::constants::math::TwoPI * radius); if (den < OneEventHitDensity(multiplicity, radius)) den = OneEventHitDensity(multiplicity, radius); @@ -301,6 +286,7 @@ float FastTracker::ProbGoodChiSqHit(float radius, float searchRadiusRPhi, float // returns number of intercepts (generic for now) int FastTracker::FastTrack(o2::track::TrackParCov inputTrack, o2::track::TrackParCov& outputTrack, const float nch) { + dNdEtaCent = nch; // set the number of charged particles per unit rapidity hits.clear(); nIntercepts = 0; nSiliconPoints = 0; @@ -310,38 +296,54 @@ int FastTracker::FastTrack(o2::track::TrackParCov inputTrack, o2::track::TrackPa const float initialRadius = std::hypot(posIni[0], posIni[1]); const float kTrackingMargin = 0.1; const int kMaxNumberOfDetectors = 20; + if (kMaxNumberOfDetectors < layers.size()) { + LOG(fatal) << "Too many layers in FastTracker, increase kMaxNumberOfDetectors"; + return -1; // too many layers + } + int firstActiveLayer = -1; // first layer that is not inert + for (size_t i = 0; i < layers.size(); ++i) { + if (!layers[i].isInert()) { + firstActiveLayer = i; + break; + } + } + if (firstActiveLayer <= 0) { + LOG(fatal) << "No active layers found in FastTracker, check layer setup"; + return -2; // no active layers + } const int xrhosteps = 100; const bool applyAngularCorrection = true; goodHitProbability.clear(); - for (int i = 0; i < kMaxNumberOfDetectors; ++i) + for (int i = 0; i < kMaxNumberOfDetectors; ++i) { goodHitProbability.push_back(-1.); - goodHitProbability[0] = 1.; + } + goodHitProbability[0] = 1.; // we use layer zero to accumulate // +-~-<*>-~-+-~-<*>-~-+-~-<*>-~-+-~-<*>-~-+-~-<*>-~-+ // Outward pass to find intercepts int firstLayerReached = -1; int lastLayerReached = -1; new (&outputTrack)(o2::track::TrackParCov)(inputTrack); - for (uint32_t il = 0; il < layers.size(); il++) { + for (size_t il = 0; il < layers.size(); il++) { // check if layer is doable - if (layers[il].getRadius() < initialRadius) + if (layers[il].getRadius() < initialRadius) { continue; // this layer should not be attempted, but go ahead + } // check if layer is reached float targetX = 1e+3; - bool ok = true; inputTrack.getXatLabR(layers[il].getRadius(), targetX, magneticField); if (targetX > 999.f) { LOGF(debug, "Failed to find intercept for layer %d at radius %.2f cm", il, layers[il].getRadius()); break; // failed to find intercept } - ok = inputTrack.propagateTo(targetX, magneticField); - if (ok && applyMSCorrection && layers[il].getRadiationLength() > 0) { + bool ok = inputTrack.propagateTo(targetX, magneticField); + if (ok && mApplyMSCorrection && layers[il].getRadiationLength() > 0) { ok = inputTrack.correctForMaterial(layers[il].getRadiationLength(), 0, applyAngularCorrection); } - if (ok && applyElossCorrection && layers[il].getDensity() > 0) { // correct in small steps + if (ok && mApplyElossCorrection && layers[il].getDensity() > 0) { // correct in small steps for (int ise = xrhosteps; ise--;) { ok = inputTrack.correctForMaterial(0, -layers[il].getDensity() / xrhosteps, applyAngularCorrection); if (!ok) @@ -361,7 +363,7 @@ int FastTracker::FastTrack(o2::track::TrackParCov inputTrack, o2::track::TrackPa break; } } - if (std::abs(inputTrack.getZ()) > layers[il].getZ() && applyZacceptance) { + if (std::abs(inputTrack.getZ()) > layers[il].getZ() && mApplyZacceptance) { break; // out of acceptance bounds } @@ -384,39 +386,19 @@ int FastTracker::FastTrack(o2::track::TrackParCov inputTrack, o2::track::TrackPa o2::track::TrackParCov inwardTrack(inputTrack); // Enlarge covariance matrix - std::array trPars = {0.}; - for (int ip = 0; ip < 5; ip++) { + std::array trPars = {0.}; + for (int ip = 0; ip < o2::track::kNParams; 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--;) + static constexpr float kLargeErr2Coord = 5 * 5; + static constexpr float kLargeErr2Dir = 0.7 * 0.7; + static constexpr float kLargeErr2PtI = 30.5 * 30.5; + std::array largeCov = {0.}; + for (int ic = o2::track::kCovMatSize; ic--;) largeCov[ic] = 0.; - largeCov[kY2] = largeCov[kZ2] = kLargeErr2Coord; - largeCov[kSnp2] = largeCov[kTgl2] = kLargeErr2Dir; - largeCov[kPtI2] = kLargeErr2PtI * trPars[kPtI] * trPars[kPtI]; + largeCov[o2::track::CovLabels::kSigY2] = largeCov[o2::track::CovLabels::kSigZ2] = kLargeErr2Coord; + largeCov[o2::track::CovLabels::kSigSnp2] = largeCov[o2::track::CovLabels::kSigTgl2] = kLargeErr2Dir; + largeCov[o2::track::CovLabels::kSigQ2Pt2] = kLargeErr2PtI * trPars[o2::track::ParLabels::kQ2Pt] * trPars[o2::track::ParLabels::kQ2Pt]; inwardTrack.setCov(largeCov); inwardTrack.checkCovariance(); @@ -434,7 +416,7 @@ int FastTracker::FastTrack(o2::track::TrackParCov inputTrack, o2::track::TrackPa continue; // failed to propagate } - if (std::abs(inputTrack.getZ()) > layers[il].getZ() && applyZacceptance) { + if (std::abs(inputTrack.getZ()) > layers[il].getZ() && mApplyZacceptance) { continue; // out of acceptance bounds but continue inwards } @@ -444,10 +426,10 @@ int FastTracker::FastTrack(o2::track::TrackParCov inputTrack, o2::track::TrackPa std::vector thisHit = {spacePoint[0], spacePoint[1], spacePoint[2]}; // towards adding cluster: move to track alpha - double alpha = inwardTrack.getAlpha(); - double xyz1[3]{ - TMath::Cos(alpha) * spacePoint[0] + TMath::Sin(alpha) * spacePoint[1], - -TMath::Sin(alpha) * spacePoint[0] + TMath::Cos(alpha) * spacePoint[1], + float alpha = inwardTrack.getAlpha(); + float xyz1[3]{ + std::cos(alpha) * spacePoint[0] + std::sin(alpha) * spacePoint[1], + -std::sin(alpha) * spacePoint[0] + std::cos(alpha) * spacePoint[1], spacePoint[2]}; if (!inwardTrack.propagateTo(xyz1[0], magneticField)) continue; @@ -462,7 +444,7 @@ int FastTracker::FastTrack(o2::track::TrackParCov inputTrack, o2::track::TrackPa inwardTrack.checkCovariance(); } - if (applyMSCorrection && layers[il].getRadiationLength() > 0) { + if (mApplyMSCorrection && layers[il].getRadiationLength() > 0) { if (!inputTrack.correctForMaterial(layers[il].getRadiationLength(), 0, applyAngularCorrection)) { return -6; } @@ -470,7 +452,7 @@ int FastTracker::FastTrack(o2::track::TrackParCov inputTrack, o2::track::TrackPa return -6; } } - if (applyElossCorrection && layers[il].getDensity() > 0) { + if (mApplyElossCorrection && layers[il].getDensity() > 0) { for (int ise = xrhosteps; ise--;) { // correct in small steps if (!inputTrack.correctForMaterial(0, layers[il].getDensity() / xrhosteps, applyAngularCorrection)) { return -7; @@ -488,9 +470,9 @@ int FastTracker::FastTrack(o2::track::TrackParCov inputTrack, o2::track::TrackPa hits.push_back(thisHit); - if (applyEffCorrection && !layers[il].isInert()) { // good hit probability calculation - double sigYCmb = o2::math_utils::sqrt(inwardTrack.getSigmaY2() + layers[il].getResolutionRPhi() * layers[il].getResolutionRPhi()); - double sigZCmb = o2::math_utils::sqrt(inwardTrack.getSigmaZ2() + layers[il].getResolutionZ() * layers[il].getResolutionZ()); + if (!layers[il].isInert()) { // good hit probability calculation + float sigYCmb = o2::math_utils::sqrt(inwardTrack.getSigmaY2() + layers[il].getResolutionRPhi() * layers[il].getResolutionRPhi()); + float sigZCmb = o2::math_utils::sqrt(inwardTrack.getSigmaZ2() + layers[il].getResolutionZ() * layers[il].getResolutionZ()); goodHitProbability[il] = ProbGoodChiSqHit(layers[il].getRadius() * 100, sigYCmb * 100, sigZCmb * 100); goodHitProbability[0] *= goodHitProbability[il]; } @@ -498,9 +480,11 @@ int FastTracker::FastTrack(o2::track::TrackParCov inputTrack, o2::track::TrackPa // backpropagate to original radius float finalX = 1e+3; - inwardTrack.getXatLabR(initialRadius, finalX, magneticField); - if (finalX > 999) + bool inPropStatus = inwardTrack.getXatLabR(initialRadius, finalX, magneticField); + if (finalX > 999) { + LOG(debug) << "Failed to find intercept for initial radius " << initialRadius << " cm, x = " << finalX << " and status " << inPropStatus << " and sn = " << inwardTrack.getSnp() << " r = " << inwardTrack.getY() * inwardTrack.getY(); return -3; // failed to find intercept + } if (!inwardTrack.propagateTo(finalX, magneticField)) { return -4; // failed to propagate @@ -511,17 +495,15 @@ int FastTracker::FastTrack(o2::track::TrackParCov inputTrack, o2::track::TrackPa return nIntercepts; // generate efficiency - if (applyEffCorrection) { - dNdEtaCent = nch; - float eff = 1.; - for (int i = 0; i < kMaxNumberOfDetectors; i++) { - float iGoodHit = goodHitProbability[i]; - if (iGoodHit <= 0) - continue; - - eff *= iGoodHit; - } + float eff = 1.; + for (int i = 0; i < kMaxNumberOfDetectors; i++) { + float iGoodHit = goodHitProbability[i]; + if (iGoodHit <= 0) + continue; + eff *= iGoodHit; + } + if (mApplyEffCorrection) { if (gRandom->Uniform() > eff) return -8; } @@ -530,11 +512,11 @@ int FastTracker::FastTrack(o2::track::TrackParCov inputTrack, o2::track::TrackPa outputTrack.checkCovariance(); // Use covariance matrix based smearing - std::array covMat = {0.}; - for (int ii = 0; ii < 15; ii++) + std::array covMat = {0.}; + for (int ii = 0; ii < o2::track::kCovMatSize; ii++) covMat[ii] = outputTrack.getCov()[ii]; TMatrixDSym m(5); - double fcovm[5][5]; + float fcovm[5][5]; for (int ii = 0, k = 0; ii < 5; ++ii) { for (int j = 0; j < ii + 1; ++j, ++k) { @@ -544,7 +526,7 @@ int FastTracker::FastTrack(o2::track::TrackParCov inputTrack, o2::track::TrackPa } // evaluate ruben's conditional, regularise - bool makePositiveDefinite = (covMatFactor > -1e-5); // apply fix + const bool makePositiveDefinite = (covMatFactor > -1e-5); // apply fix bool rubenConditional = false; for (int ii = 0; ii < 5; ii++) { for (int jj = 0; jj < 5; jj++) { @@ -571,7 +553,7 @@ int FastTracker::FastTrack(o2::track::TrackParCov inputTrack, o2::track::TrackPa } if (negEigVal && rubenConditional && makePositiveDefinite) { - if (verboseLevel > 0) { + if (mVerboseLevel > 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; @@ -585,9 +567,9 @@ int FastTracker::FastTrack(o2::track::TrackParCov inputTrack, o2::track::TrackPa covMatOK++; // transform parameter vector and smear - double params_[5]; + float params_[5]; for (int ii = 0; ii < 5; ++ii) { - double val = 0.; + float val = 0.; for (int j = 0; j < 5; ++j) val += eigVec[j][ii] * outputTrack.getParam(j); // smear parameters according to eigenvalues @@ -598,7 +580,7 @@ int FastTracker::FastTrack(o2::track::TrackParCov inputTrack, o2::track::TrackPa eigVec.Invert(); // transform back params vector for (int ii = 0; ii < 5; ++ii) { - double val = 0.; + float val = 0.; for (int j = 0; j < 5; ++j) val += eigVec[j][ii] * params_[j]; outputTrack.setParam(val, ii); diff --git a/ALICE3/Core/FastTracker.h b/ALICE3/Core/FastTracker.h index 702caa9e84b..a0dba5d7ec5 100644 --- a/ALICE3/Core/FastTracker.h +++ b/ALICE3/Core/FastTracker.h @@ -12,12 +12,15 @@ #ifndef ALICE3_CORE_FASTTRACKER_H_ #define ALICE3_CORE_FASTTRACKER_H_ -#include // not a system header but megalinter thinks so -#include -#include #include "DetLayer.h" + #include "ReconstructionDataFormats/Track.h" +#include // not a system header but megalinter thinks so + +#include +#include + namespace o2 { namespace fastsim @@ -32,13 +35,16 @@ class FastTracker { public: // Constructor/destructor - FastTracker(); + FastTracker() = default; + // Destructor virtual ~FastTracker() {} // Layer and layer configuration void AddLayer(TString name, float r, float z, float x0, float xrho, float resRPhi = 0.0f, float resZ = 0.0f, float eff = 0.0f, int type = 0); DetLayer GetLayer(const int layer, bool ignoreBarrelLayers = true) const; - int GetLayerIndex(const std::string name) const; + int GetLayerIndex(const std::string& name) const; + size_t GetNLayers() const { return layers.size(); } + bool IsLayerInert(const int layer) const { return layers[layer].isInert(); } void SetRadiationLength(const std::string layerName, float x0) { layers[GetLayerIndex(layerName)].setRadiationLength(x0); } void SetRadius(const std::string layerName, float r) { layers[GetLayerIndex(layerName)].setRadius(r); } void SetResolutionRPhi(const std::string layerName, float resRPhi) { layers[GetLayerIndex(layerName)].setResolutionRPhi(resRPhi); } @@ -80,20 +86,24 @@ class FastTracker void SetIntegrationTime(float t) { integrationTime = t; } void SetMaxRadiusOfSlowDetectors(float r) { maxRadiusSlowDet = r; } void SetAvgRapidity(float y) { avgRapidity = y; } - void SetdNdEtaCent(float d) { dNdEtaCent = d; } + void SetdNdEtaCent(int d) { dNdEtaCent = d; } void SetLhcUPCscale(float s) { lhcUPCScale = s; } void SetBField(float b) { magneticField = b; } void SetMinRadTrack(float r) { fMinRadTrack = r; } void SetMagneticField(float b) { magneticField = b; } - void SetApplyZacceptance(bool b) { applyZacceptance = b; } - void SetApplyMSCorrection(bool b) { applyMSCorrection = b; } - void SetApplyElossCorrection(bool b) { applyElossCorrection = b; } + void SetApplyZacceptance(bool b) { mApplyZacceptance = b; } + void SetApplyMSCorrection(bool b) { mApplyMSCorrection = b; } + void SetApplyElossCorrection(bool b) { mApplyElossCorrection = b; } + void SetApplyEffCorrection(bool b) { mApplyEffCorrection = b; } // Getters for the last track int GetNIntercepts() const { return nIntercepts; } int GetNSiliconPoints() const { return nSiliconPoints; } int GetNGasPoints() const { return nGasPoints; } - float GetGoodHitProb(int layer) const { return goodHitProbability[layer]; } + float GetGoodHitProb(int layer) const + { + return (layer >= 0 && static_cast(layer) < goodHitProbability.size()) ? goodHitProbability[layer] : 0.0f; + } std::size_t GetNHits() const { return hits.size(); } float GetHitX(const int i) const { return hits[i][0]; } float GetHitY(const int i) const { return hits[i][1]; } @@ -106,34 +116,35 @@ class FastTracker std::vector layers; std::vector> hits; // bookkeep last added hits - // operational - bool applyZacceptance; // check z acceptance or not - bool applyMSCorrection; // Apply correction for multiple scattering - bool applyElossCorrection; // Apply correction for eloss (requires MS correction) - bool applyEffCorrection; // Apply correction for hit efficiency - int verboseLevel; // 0: not verbose, >0 more verbose - int crossSectionMinB; - int dNdEtaCent; - int dNdEtaMinB; - float integrationTime; - float magneticField; // in kiloGauss (5 = 0.5T, etc) - float covMatFactor; // covmat off-diagonal factor to use for covmat fix (negative: no factor) - float sigmaD; - float luminosity; - float otherBackground; - float maxRadiusSlowDet; - float avgRapidity; - float lhcUPCScale; - float upcBackgroundMultiplier; - float fMinRadTrack = 132.; - - uint64_t covMatOK; // cov mat has negative eigenvals - uint64_t covMatNotOK; // cov mat has negative eigenvals - - // last track information - int nIntercepts; // found in first outward propagation - int nSiliconPoints; // silicon-based space points added to track - int nGasPoints; // tpc-based space points added to track + /// configuration parameters + bool mApplyZacceptance = false; /// check z acceptance or not + bool mApplyMSCorrection = true; /// Apply correction for multiple scattering + bool mApplyElossCorrection = true; /// Apply correction for eloss (requires MS correction) + bool mApplyEffCorrection = true; /// Apply correction for hit efficiency + int mVerboseLevel = 0; /// 0: not verbose, >0 more verbose + const float mCrossSectionMinB = 8; /// Minimum bias Cross section for event under study (PbPb MinBias ~ 8 Barns) + int dNdEtaCent = 2200; /// dN/deta e.g. at centrality 0-5% (for 5 TeV PbPb) + int dNdEtaMinB = 1; /// dN/deta for minimum bias events + float integrationTime = 0.02f; /// Integration time in ms + float magneticField = 20.f; /// Magnetic field in kiloGauss (5 = 0.5T, 20 = 2T, etc) + float covMatFactor = 0.99f; /// covmat off-diagonal factor to use for covmat fix (negative: no factor) + float sigmaD = 6.0f; /// sigma for the detector resolution in cm + float luminosity = 1.e27f; /// luminosity in cm^-2 s^-1 (e.g. 1.e27 for PbPb at 5 TeV) + float otherBackground = 0.0f; /// background from other sources, e.g. pileup, in [0, 1] + float maxRadiusSlowDet = 10.f; /// maximum radius of slow detectors in cm + float avgRapidity = 0.45f; /// average rapidity for hit density calculation + float lhcUPCScale = 1.0f; /// scale factor for LHC UPC events + float upcBackgroundMultiplier = 1.0f; /// multiplier for UPC background + float fMinRadTrack = 132.f; /// minimum radius for track propagation in cm + + /// counters for covariance matrix statuses + uint64_t covMatOK = 0; /// cov mat has positive eigenvals + uint64_t covMatNotOK = 0; /// cov mat has negative eigenvals + + /// last track information + int nIntercepts = 0; /// found in first outward propagation + int nSiliconPoints = 0; /// silicon-based space points added to track + int nGasPoints = 0; /// tpc-based space points added to track std::vector goodHitProbability; ClassDef(FastTracker, 1);