diff --git a/DataFormats/Detectors/TRD/include/DataFormatsTRD/CalGain.h b/DataFormats/Detectors/TRD/include/DataFormatsTRD/CalGain.h index f90101e7a4f21..408198665df6a 100644 --- a/DataFormats/Detectors/TRD/include/DataFormatsTRD/CalGain.h +++ b/DataFormats/Detectors/TRD/include/DataFormatsTRD/CalGain.h @@ -33,12 +33,42 @@ class CalGain void setMPVdEdx(int iDet, float mpv) { mMPVdEdx[iDet] = mpv; } - float getMPVdEdx(int iDet) const { return mMPVdEdx[iDet]; } + float getMPVdEdx(int iDet, bool defaultAvg = false) const { + // if defaultAvg = false, we take the value stored whatever it is + // if defaultAvg = true and we have default value or bad value stored, we take the average on all chambers instead + if (!defaultAvg || isGoodGain(iDet)) return mMPVdEdx[iDet]; + else return getAverageGain(); + } + + + float getAverageGain() const { + float averageGain = 0.; + int ngood = 0; + + for (int iDet = 0; iDet < constants::MAXCHAMBER; iDet++) { + if (isGoodGain(iDet)) { + // The chamber has correct calibration + ngood ++; + averageGain += mMPVdEdx[iDet]; + } + } + if (ngood == 0) { + // we should make sure it never happens + return constants::MPVDEDXDEFAULT; + } + averageGain /= ngood; + return averageGain; + } + + bool isGoodGain(int iDet) const { + if (TMath::Abs(mMPVdEdx[iDet] - constants::MPVDEDXDEFAULT) > 1e-6) return true; + else return false; + } private: std::array mMPVdEdx{}; ///< Most probable value of dEdx distribution per TRD chamber - ClassDefNV(CalGain, 1); + ClassDefNV(CalGain, 2); }; } // namespace trd diff --git a/DataFormats/Detectors/TRD/include/DataFormatsTRD/CalVdriftExB.h b/DataFormats/Detectors/TRD/include/DataFormatsTRD/CalVdriftExB.h index bad9dcfef4e37..aa47a75a65e7a 100644 --- a/DataFormats/Detectors/TRD/include/DataFormatsTRD/CalVdriftExB.h +++ b/DataFormats/Detectors/TRD/include/DataFormatsTRD/CalVdriftExB.h @@ -33,15 +33,83 @@ class CalVdriftExB void setVdrift(int iDet, float vd) { mVdrift[iDet] = vd; } void setExB(int iDet, float exb) { mExB[iDet] = exb; } + + float getVdrift(int iDet, bool defaultAvg = false) const { + // if defaultAvg = false, we take the value stored whatever it is + // if defaultAvg = true and we have default value or bad value stored, we take the average on all chambers instead + if (!defaultAvg || (isGoodExB(iDet) && isGoodVdrift(iDet))) return mVdrift[iDet]; + else return getAverageVdrift(); + } + float getExB(int iDet, bool defaultAvg = false) const { + if (!defaultAvg || (isGoodExB(iDet) && isGoodVdrift(iDet))) return mExB[iDet]; + else return getAverageExB(); + } + + float getAverageVdrift() const { + float averageVdrift = 0.; + int ngood = 0; + + for (int iDet = 0; iDet < constants::MAXCHAMBER; iDet++) { + if (isGoodExB(iDet) && isGoodVdrift(iDet)) { + // Both values need to be correct to declare a chamber as well calibrated + ngood ++; + averageVdrift += mVdrift[iDet]; + } + } + if (ngood == 0) { + // we should make sure it never happens + return constants::VDRIFTDEFAULT; + } + averageVdrift /= ngood; + return averageVdrift; + } - float getVdrift(int iDet) const { return mVdrift[iDet]; } - float getExB(int iDet) const { return mExB[iDet]; } + float getAverageExB() const { + float averageExB = 0.; + int ngood = 0; + + for (int iDet = 0; iDet < constants::MAXCHAMBER; iDet++) { + if (isGoodExB(iDet) && isGoodVdrift(iDet)) { + // Both values need to be correct to declare a chamber as well calibrated + ngood ++; + averageExB += mExB[iDet]; + } + } + if (ngood == 0) { + // we should make sure it never happens + return constants::EXBDEFAULT; + } + averageExB /= ngood; + return averageExB; + } + + bool isGoodExB(int iDet) const { + // check if value is well calibrated or not + // default calibration if not enough entries + // close to boundaries indicate a failed fit + if (TMath::Abs(mExB[iDet] - constants::EXBDEFAULT) > 1e-6 && + TMath::Abs(mExB[iDet] - constants::EXBMIN) > 0.01 && + TMath::Abs(mExB[iDet] - constants::EXBMAX) > 0.01) + return true; + else return false; + } + + bool isGoodVdrift(int iDet) const { + // check if value is well calibrated or not + // default calibration if not enough entries + // close to boundaries indicate a failed fit + if (TMath::Abs(mVdrift[iDet] - constants::VDRIFTDEFAULT) > 1e-6 && + TMath::Abs(mVdrift[iDet] - constants::VDRIFTMIN) > 0.1 && + TMath::Abs(mVdrift[iDet] - constants::VDRIFTMAX) > 0.1) + return true; + else return false; + } private: std::array mVdrift{}; ///< calibrated drift velocity per TRD chamber std::array mExB{}; ///< calibrated Lorentz angle per TRD chamber - ClassDefNV(CalVdriftExB, 1); + ClassDefNV(CalVdriftExB, 2); }; } // namespace trd diff --git a/DataFormats/Detectors/TRD/include/DataFormatsTRD/Constants.h b/DataFormats/Detectors/TRD/include/DataFormatsTRD/Constants.h index 7a650cf3699cf..b2c827710356d 100644 --- a/DataFormats/Detectors/TRD/include/DataFormatsTRD/Constants.h +++ b/DataFormats/Detectors/TRD/include/DataFormatsTRD/Constants.h @@ -75,7 +75,11 @@ constexpr int TIMEBINS = 30; ///< the number of time bins constexpr float MAXIMPACTANGLE = 25.f; ///< the maximum impact angle for tracks relative to the TRD detector plane to be considered for vDrift and ExB calibration constexpr int NBINSANGLEDIFF = 25; ///< the number of bins for the track angle used for the vDrift and ExB calibration based on the tracking constexpr double VDRIFTDEFAULT = 1.546; ///< default value for vDrift +constexpr double VDRIFTMIN = 0.4; ///< min value for vDrift +constexpr double VDRIFTMAX = 2.0; ///< max value for vDrift constexpr double EXBDEFAULT = 0.0; ///< default value for LorentzAngle +constexpr double EXBMIN = -0.4; ///< min value for LorentzAngle +constexpr double EXBMAX = 0.4; ///< max value for LorentzAngle constexpr int NBINSGAINCALIB = 320; ///< number of bins in the charge (Q0+Q1+Q2) histogram for gain calibration constexpr float MPVDEDXDEFAULT = 42.; ///< default Most Probable Value of TRD dEdx constexpr float T0DEFAULT = 1.2; ///< default value for t0 diff --git a/Detectors/Align/Workflow/src/BarrelAlignmentSpec.cxx b/Detectors/Align/Workflow/src/BarrelAlignmentSpec.cxx index d4ab53c8181ce..c4367ef038cbc 100644 --- a/Detectors/Align/Workflow/src/BarrelAlignmentSpec.cxx +++ b/Detectors/Align/Workflow/src/BarrelAlignmentSpec.cxx @@ -307,7 +307,12 @@ void BarrelAlignmentSpec::finaliseCCDB(o2::framework::ConcreteDataMatcher& match } if (matcher == ConcreteDataMatcher("TRD", "CALVDRIFTEXB", 0)) { LOG(info) << "CalVdriftExB object has been updated"; - mTRDTransformer->setCalVdriftExB((const o2::trd::CalVdriftExB*)obj); + for (int iDet = 0; iDet < o2::trd::constants::MAXCHAMBER; iDet++) { + // set to average value if the calibration is not correct + mTRDTransformer->setVdrift(iDet, ((const o2::trd::CalVdriftExB*)obj)->getVdrift(iDet, true)); + mTRDTransformer->setExB(iDet, ((const o2::trd::CalVdriftExB*)obj)->getExB(iDet, true)); + } + //mTRDTransformer->setCalVdriftExB((const o2::trd::CalVdriftExB*)obj); return; } if (mTPCVDriftHelper.accountCCDBInputs(matcher, obj)) { diff --git a/Detectors/TRD/base/include/TRDBase/TrackletTransformer.h b/Detectors/TRD/base/include/TRDBase/TrackletTransformer.h index ca9e71d392231..1defe41ec9f58 100644 --- a/Detectors/TRD/base/include/TRDBase/TrackletTransformer.h +++ b/Detectors/TRD/base/include/TRDBase/TrackletTransformer.h @@ -30,10 +30,13 @@ class TrackletTransformer void init(); - void setCalVdriftExB(const CalVdriftExB* cal) { mCalVdriftExB = cal; }; + void setVdrift(int iDet, float vd) {mVdrift[iDet] = vd;} + void setExB(int iDet, float exb) {mExB[iDet] = exb;} void setApplyXOR() { mApplyXOR = true; } void setApplyShift(bool f) { mApplyShift = f; } bool isShiftApplied() const { return mApplyShift; } + + float calculateZ(int padrow, const PadPlane* padPlane) const; @@ -54,7 +57,8 @@ class TrackletTransformer float mXAnode; - const CalVdriftExB* mCalVdriftExB{nullptr}; + std::array mVdrift{}; + std::array mExB{}; }; } // namespace trd diff --git a/Detectors/TRD/base/src/TrackletTransformer.cxx b/Detectors/TRD/base/src/TrackletTransformer.cxx index 58938cfd99161..370426d7c7201 100644 --- a/Detectors/TRD/base/src/TrackletTransformer.cxx +++ b/Detectors/TRD/base/src/TrackletTransformer.cxx @@ -40,8 +40,10 @@ float TrackletTransformer::calculateDy(int detector, int slope, const PadPlane* { double padWidth = padPlane->getWidthIPad(); - float vDrift = mCalVdriftExB->getVdrift(detector); - float exb = mCalVdriftExB->getExB(detector); + //float vDrift = mCalVdriftExB->getVdrift(detector, true); + //float exb = mCalVdriftExB->getExB(detector, true); + float vDrift = mVdrift[detector]; + float exb = mExB[detector]; // dy = slope * nTimeBins * padWidth * GRANULARITYTRKLSLOPE; // nTimeBins should be number of timebins in drift region. 1 timebin is 100 nanosecond @@ -50,11 +52,12 @@ float TrackletTransformer::calculateDy(int detector, int slope, const PadPlane* // NOTE: check what drift height is used in calibration code to ensure consistency // NOTE: check sign convention of Lorentz angle // NOTE: confirm the direction in which vDrift is measured/determined. Is it in x or in direction of drift? - double lorentzCorrection = TMath::Tan(exb) * mXAnode; + // The Lorentz correction have to be applied both at the point of entrance and at the end of the drift region + double lorentzCorrection = TMath::Tan(exb) * mGeo->cdrHght(); // assuming angle in Bailhache, fig. 4.17 would be positive in our calibration code double calibratedDy = rawDy - lorentzCorrection; - + return calibratedDy; } @@ -96,7 +99,7 @@ CalibratedTracklet TrackletTransformer::transformTracklet(Tracklet64 tracklet, b position = tracklet.getPositionBinSigned(); slope = tracklet.getSlopeBinSigned(); } - + // calculate raw local chamber space point const auto padPlane = mGeo->getPadPlane(detector); @@ -124,7 +127,7 @@ CalibratedTracklet TrackletTransformer::transformTracklet(Tracklet64 tracklet, b double TrackletTransformer::getTimebin(int detector, double x) const { // calculate timebin from x position within chamber - float vDrift = mCalVdriftExB->getVdrift(detector); + float vDrift = mVdrift[detector]; double t0 = 4.0; // time (in timebins) of start of drift region double timebin; diff --git a/Detectors/TRD/calibration/README.md b/Detectors/TRD/calibration/README.md index bdbfc9e709800..f769fa99b5778 100644 --- a/Detectors/TRD/calibration/README.md +++ b/Detectors/TRD/calibration/README.md @@ -27,15 +27,15 @@ For 'o2-calibration-trd-workflow --vDriftAndExB' there are also the following ke *Hint: You can get information on the meaning of the parameters by running `o2-calibration-trd-workflow --vDriftAndExB -b --help full`* -If you want to run the calibration from a local file with residuals, trdangreshistos.root, you can run: +If you want to run the calibration from a local file with residuals, trdcaliboutput.root, you can run: - o2-calibration-trd-workflow --vDriftAndExB -b --enable-root-input --calib-vdexb-calibration '--tf-per-slot 1' --configKeyValues "TRDCalibParams.minEntriesChamber=100;TRDCalibParams.minEntriesTotal=50000" + o2-calibration-trd-workflow --vDriftAndExB -b --enable-root-input --calib-vdexb-calibration '--tf-per-slot 1' --configKeyValues "TRDCalibParams.minEntriesChamber=100;TRDCalibParams.minEntriesTotal=50000" --trd-calib-infile trdcaliboutput.root Additionally it is possible to perform the calibrations fit manually per chamber if you have TPC-TRD or ITS-TPC-TRD tracks, you can run: o2-trd-global-tracking -b --enable-trackbased-calib -This produces `trdangreshistos.root` which holds the residuals of the angles and differences. +This produces `trdcaliboutput.root` which holds the residuals of the angles and differences. Then run the macro `Detectors/TRD/calibration/macros/manualCalibFit.C`. This produces a file of similar name with the fitted data and prints out the fit results. This is equivalent to running: diff --git a/Detectors/TRD/calibration/include/TRDCalibration/CalibrationParams.h b/Detectors/TRD/calibration/include/TRDCalibration/CalibrationParams.h index 677673a7f85f3..31f7f56850419 100644 --- a/Detectors/TRD/calibration/include/TRDCalibration/CalibrationParams.h +++ b/Detectors/TRD/calibration/include/TRDCalibration/CalibrationParams.h @@ -25,12 +25,13 @@ namespace trd /// VDrift and ExB calibration parameters. struct TRDCalibParams : public o2::conf::ConfigurableParamHelper { unsigned int nTrackletsMin = 5; ///< minimum amount of tracklets + unsigned int nTrackletsMinLoose = 4; ///< minimum amount of tracklets if two layers with a large lever arm both have a hit unsigned int chi2RedMax = 6; ///< maximum reduced chi2 acceptable for track quality - size_t minEntriesChamber = 75; ///< minimum number of entries per chamber to fit single time slot - size_t minEntriesTotal = 40'500; ///< minimum total required for meaningful fits + size_t minEntriesChamber = 200; ///< minimum number of entries per chamber to fit single time slot + size_t minEntriesTotal = 400'000; ///< minimum total required for meaningful fits // For gain calibration - unsigned int nTrackletsMinGainCalib = 5; + unsigned int nTrackletsMinGainCalib = 3; size_t minEntriesChamberGainCalib = 500; ///< minimum number of entries per chamber to fit single time slot size_t minEntriesTotalGainCalib = 1'000'000; ///< minimum total required for meaningful fits // Cuts for selecting clean pion candidates for gain calibration diff --git a/Detectors/TRD/calibration/include/TRDCalibration/CalibratorVdExB.h b/Detectors/TRD/calibration/include/TRDCalibration/CalibratorVdExB.h index 7d55850af9fd4..50cc3930baa96 100644 --- a/Detectors/TRD/calibration/include/TRDCalibration/CalibratorVdExB.h +++ b/Detectors/TRD/calibration/include/TRDCalibration/CalibratorVdExB.h @@ -88,8 +88,8 @@ class CalibratorVdExB final : public o2::calibration::TimeSlotCalibration mOutFile{nullptr}; ///< output file std::unique_ptr mOutTree{nullptr}; ///< output tree diff --git a/Detectors/TRD/calibration/macros/manualCalibFit.C b/Detectors/TRD/calibration/macros/manualCalibFit.C index d31744a2e727c..c9eb5ec0d66a6 100644 --- a/Detectors/TRD/calibration/macros/manualCalibFit.C +++ b/Detectors/TRD/calibration/macros/manualCalibFit.C @@ -30,18 +30,19 @@ // O2 header #include +#include "DetectorsBase/Propagator.h" #endif // This root macro reads in 'trdangreshistos.root' and // performs the calibration fits manually as in CalibratorVdExB.cxx // This can be used for checking if the calibration fits make sense. -void manualCalibFit() +void manualCalibFit(int runNumber = 563335, bool usePreCorrFromCCDB = false) { //---------------------------------------------------- // TTree and File //---------------------------------------------------- - std::unique_ptr inFilePtr(TFile::Open("trdangreshistos.root")); + std::unique_ptr inFilePtr(TFile::Open("trdcaliboutput.root")); if (inFilePtr == nullptr) { printf("Input File could not be read!\n'"); return; @@ -59,19 +60,49 @@ void manualCalibFit() mNEntriesPerBinSum.fill(0); tree->SetBranchAddress("mHistogramEntries[13500]", &mHistogramEntries); tree->SetBranchAddress("mNEntriesPerBin[13500]", &mNEntriesPerBin); + + // use precorr values from ccdb + // necessary when the angular residuals were calculated already using ccdb calibration (e.g. in a local run) + + o2::trd::CalVdriftExB* calObject; + if (usePreCorrFromCCDB) { + auto& ccdbmgr = o2::ccdb::BasicCCDBManager::instance(); + + o2::ccdb::CcdbApi ccdb; + ccdb.init("http://alice-ccdb.cern.ch"); + auto runDuration = ccdbmgr.getRunDuration(runNumber); + + std::map metadata; + std::map headers; + + calObject = ccdb.retrieveFromTFileAny("TRD/Calib/CalVdriftExB", metadata, runDuration.first + 60000, &headers, "", "", "1689478811721"); + } + //---------------------------------------------------- // Configure Fitter //---------------------------------------------------- o2::trd::FitFunctor mFitFunctor; std::array, 540> profiles; ///< profile histograms for each TRD chamber + int counter = 0; for (int iDet = 0; iDet < 540; ++iDet) { mFitFunctor.profiles[iDet] = std::make_unique(Form("profAngleDiff_%i", iDet), Form("profAngleDiff_%i", iDet), 25, -25.f, 25.f); + if (usePreCorrFromCCDB) { + std::cout<getVdrift(iDet)<<" "<getExB(iDet)<<" "; + if (calObject->isGoodExB(iDet)) counter++; + if (iDet%6==5)std::cout<getVdrift(iDet, true); + mFitFunctor.laPreCorr[iDet] = calObject->getExB(iDet, true); + } } + std::cout << counter << " good entries in the CCDB " << std::endl; + mFitFunctor.mAnodePlane = 3.35; // don't really care as long as it's not zero, this parameter could be removed mFitFunctor.lowerBoundAngleFit = 80 * TMath::DegToRad(); mFitFunctor.upperBoundAngleFit = 100 * TMath::DegToRad(); - mFitFunctor.vdPreCorr.fill(1.546); - mFitFunctor.laPreCorr.fill(0.0); + if (!usePreCorrFromCCDB) { + mFitFunctor.vdPreCorr.fill(1.546); + mFitFunctor.laPreCorr.fill(0.0); + } //---------------------------------------------------- // Loop @@ -88,15 +119,19 @@ void manualCalibFit() //---------------------------------------------------- // Fill profiles //---------------------------------------------------- + int nEntriesDetTotal[540] = {}; for (int iDet = 0; iDet < 540; ++iDet) { for (int iBin = 0; iBin < 25; ++iBin) { auto angleDiffSum = mHistogramEntriesSum[iDet * 25 + iBin]; auto nEntries = mNEntriesPerBinSum[iDet * 25 + iBin]; + nEntriesDetTotal[iDet] += nEntries; if (nEntries > 0) { // skip entries which have no entries; ? // add to the respective profile for fitting later on mFitFunctor.profiles[iDet]->Fill(2 * iBin - 25.f, angleDiffSum / nEntries, nEntries); + if (iDet == 207) std::cout<GetBinEntries(iBin+1)<<" "<GetBinEffectiveEntries(iBin+1)<<" "<GetBinError(iBin+1)< laFitResults{}; std::array vdFitResults{}; + + TH1F* hVd = new TH1F("hVd", "v drift", 150, 0.5, 2.); + TH1F* hLa = new TH1F("hLa", "lorentz angle", 200, -25., 25.); + o2::trd::CalVdriftExB* calObjectOut = new o2::trd::CalVdriftExB(); + for (int iDet = 0; iDet < 540; ++iDet) { + if (nEntriesDetTotal[iDet] < 75) continue; mFitFunctor.currDet = iDet; ROOT::Fit::Fitter fitter; double paramsStart[2]; - paramsStart[0] = 0. * TMath::DegToRad(); + paramsStart[0] = 0.; paramsStart[1] = 1.; fitter.SetFCN(2, mFitFunctor, paramsStart); fitter.Config().ParSettings(0).SetLimits(-0.7, 0.7); fitter.Config().ParSettings(0).SetStepSize(.01); - fitter.Config().ParSettings(1).SetLimits(0., 3.); + fitter.Config().ParSettings(1).SetLimits(0.01, 3.); fitter.Config().ParSettings(1).SetStepSize(.01); ROOT::Math::MinimizerOptions opt; opt.SetMinimizerType("Minuit2"); @@ -127,14 +168,28 @@ void manualCalibFit() auto fitResult = fitter.Result(); laFitResults[iDet] = fitResult.Parameter(0); vdFitResults[iDet] = fitResult.Parameter(1); - printf("Det %d: la=%f\tvd=%f\n", iDet, laFitResults[iDet] * TMath::RadToDeg(), vdFitResults[iDet]); + if (fitResult.MinFcnValue() > 0.03) continue; + printf("Det %d: la=%.3f \tvd=%.3f \t100*minValue=%f \tentries=%d\n", iDet, laFitResults[iDet] * TMath::RadToDeg(), vdFitResults[iDet], 100*fitResult.MinFcnValue(), nEntriesDetTotal[iDet]); + hVd->Fill(vdFitResults[iDet]); + hLa->Fill(laFitResults[iDet]* TMath::RadToDeg()); + calObjectOut->setVdrift(iDet, vdFitResults[iDet]); + calObjectOut->setExB(iDet, laFitResults[iDet]); } printf("-------- Finished fits\n"); + + std::cout<<"number of chambers with enough entries: "<GetEntries()<GetMean()<<" sigma: "<GetStdDev()<GetMean()<<" sigma: "<GetStdDev()< outFilePtr(TFile::Open("manualCalibFit.root", "RECREATE")); + hVd->Write(); + hLa->Write(); + outFilePtr->WriteObjectAny(calObjectOut, "o2::trd::CalVdriftExB", "calObject"); for (auto& p : mFitFunctor.profiles) p->Write(); + } diff --git a/Detectors/TRD/calibration/src/CalibratorVdExB.cxx b/Detectors/TRD/calibration/src/CalibratorVdExB.cxx index 64a8664640e41..68f2998dd26a2 100644 --- a/Detectors/TRD/calibration/src/CalibratorVdExB.cxx +++ b/Detectors/TRD/calibration/src/CalibratorVdExB.cxx @@ -24,6 +24,7 @@ #include "CCDB/BasicCCDBManager.h" #include "CommonUtils/NameConf.h" #include "CommonUtils/MemFileHelper.h" +//#include "DetectorsBase/Propagator.h" #include #include @@ -105,17 +106,30 @@ void CalibratorVdExB::initProcessing() return; } - mFitFunctor.lowerBoundAngleFit = 80 * TMath::DegToRad(); - mFitFunctor.upperBoundAngleFit = 100 * TMath::DegToRad(); + // fit is done in region where ion tails are small, close to lorentz angle + // we want an approximate value of the lorentz angle in order to define better fit boundaries + // TODO: find a way to obtain the magnetic field even in standalone calibration + //float bz = o2::base::Propagator::Instance()->getNominalBz(); + // default angle with zero field is slightly shifted + float lorentzAngleAvg = -1.f; + /*if (TMath::Abs(bz - 2) < 0.1f) { lorentzAngleAvg = 2.f;} + if (TMath::Abs(bz + 2) < 0.1f) { lorentzAngleAvg = -4.f;} + if (TMath::Abs(bz - 5) < 0.1f) { lorentzAngleAvg = 7.f;} + if (TMath::Abs(bz + 5) < 0.1f) { lorentzAngleAvg = -9.5f;} + + LOGP(info, "b field: {} lorentz angle start: {}", bz, lorentzAngleAvg);*/ + + mFitFunctor.lowerBoundAngleFit = (80 + lorentzAngleAvg) * TMath::DegToRad(); + mFitFunctor.upperBoundAngleFit = (100 + lorentzAngleAvg) * TMath::DegToRad(); mFitFunctor.mAnodePlane = GeometryBase::camHght() / (2.f * 100.f); for (int iDet = 0; iDet < MAXCHAMBER; ++iDet) { mFitFunctor.profiles[iDet] = std::make_unique(Form("profAngleDiff_%i", iDet), Form("profAngleDiff_%i", iDet), NBINSANGLEDIFF, -MAXIMPACTANGLE, MAXIMPACTANGLE); } mFitter.SetFCN(2, mFitFunctor, mParamsStart); - mFitter.Config().ParSettings(ParamIndex::LA).SetLimits(-0.7, 0.7); + mFitter.Config().ParSettings(ParamIndex::LA).SetLimits(constants::EXBMIN, constants::EXBMAX); mFitter.Config().ParSettings(ParamIndex::LA).SetStepSize(.01); - mFitter.Config().ParSettings(ParamIndex::VD).SetLimits(0.01, 3.); + mFitter.Config().ParSettings(ParamIndex::VD).SetLimits(constants::VDRIFTMIN, constants::VDRIFTMAX); mFitter.Config().ParSettings(ParamIndex::VD).SetStepSize(.01); ROOT::Math::MinimizerOptions opt; opt.SetMinimizerType("Minuit2"); @@ -184,17 +198,30 @@ void CalibratorVdExB::finalizeSlot(Slot& slot) } // Check if we have the minimum amount of entries if (sumEntries < mMinEntriesChamber) { - LOGF(debug, "Chamber %d did not reach minimum amount of entries for refit", iDet); + LOGF(debug, "Chamber %d did not reach minimum amount of entries for refit: %d", iDet, sumEntries); continue; } + float laPreCorrTemp = mFitFunctor.laPreCorr[iDet]; + float vdPreCorrTemp = mFitFunctor.vdPreCorr[iDet]; + // Here we start from uncalibrated values, otherwise online calibration does not work properly + mFitFunctor.laPreCorr[iDet] = EXBDEFAULT; + mFitFunctor.vdPreCorr[iDet] = VDRIFTDEFAULT; + // Reset Start Parameter mParamsStart[ParamIndex::LA] = 0.0; mParamsStart[ParamIndex::VD] = 1.0; mFitter.FitFCN(); auto fitResult = mFitter.Result(); + if (fitResult.MinFcnValue() > 0.03) { + LOGF(debug, "Chamber %d fit did not converge properly, minimization value too high: %f", iDet, fitResult.MinFcnValue()); + // The fit did not work properly, so we keep previous values + mFitFunctor.laPreCorr[iDet] = laPreCorrTemp; + mFitFunctor.vdPreCorr[iDet] = vdPreCorrTemp; + continue; + } laFitResults[iDet] = fitResult.Parameter(ParamIndex::LA); vdFitResults[iDet] = fitResult.Parameter(ParamIndex::VD); - LOGF(debug, "Fit result for chamber %i: vd=%f, la=%f", iDet, vdFitResults[iDet], laFitResults[iDet] * TMath::RadToDeg()); + LOGF(debug, "Fit result for chamber %i: vd=%f, la=%f, minimizer value=%f", iDet, vdFitResults[iDet], laFitResults[iDet] * TMath::RadToDeg(), fitResult.MinFcnValue()); // Update fit values for next fit mFitFunctor.laPreCorr[iDet] = laFitResults[iDet]; mFitFunctor.vdPreCorr[iDet] = vdFitResults[iDet]; @@ -222,7 +249,7 @@ void CalibratorVdExB::finalizeSlot(Slot& slot) auto flName = o2::ccdb::CcdbApi::generateFileName(clName); std::map metadata; // TODO: do we want to store any meta data? long startValidity = slot.getStartTimeMS() - 10 * o2::ccdb::CcdbObjectInfo::SECOND; - mInfoVector.emplace_back("TRD/Calib/CalVdriftExB", clName, flName, metadata, startValidity, startValidity + o2::ccdb::CcdbObjectInfo::HOUR); + mInfoVector.emplace_back("TRD/Calib/CalVdriftExB", clName, flName, metadata, startValidity, startValidity + 3 * o2::ccdb::CcdbObjectInfo::DAY); mObjectVector.push_back(calObject); } diff --git a/Detectors/TRD/calibration/src/TrackBasedCalib.cxx b/Detectors/TRD/calibration/src/TrackBasedCalib.cxx index 011a888a47618..e794e14873d1f 100644 --- a/Detectors/TRD/calibration/src/TrackBasedCalib.cxx +++ b/Detectors/TRD/calibration/src/TrackBasedCalib.cxx @@ -187,7 +187,10 @@ int TrackBasedCalib::doTrdOnlyTrackFits(gsl::span& tracks) for (const auto& trkIn : tracks) { if (trkIn.getNtracklets() < params.nTrackletsMin) { // with less than 3 tracklets the TRD-only refit not meaningful - continue; + if (trkIn.getNtracklets() < params.nTrackletsMinLoose || !((trkIn.getTrackletIndex(0) >= 0 && (trkIn.getTrackletIndex(NLAYER-1) >= 0 || trkIn.getTrackletIndex(NLAYER-2) >= 0))) || (trkIn.getTrackletIndex(1) >= 0 && trkIn.getTrackletIndex(NLAYER-1) >= 0)) { + // we check if we have enough lever arm, i.e. (first and last) or (second and last) or (first and before last) are present + continue; + } } auto trkWork = trkIn; // input is const, so we need to create a copy bool trackFailed = false; @@ -259,9 +262,20 @@ int TrackBasedCalib::doTrdOnlyTrackFits(gsl::span& tracks) } float trkAngle = o2::math_utils::asin(trkWork.getSnp()) * TMath::RadToDeg(); - float trkltAngle = o2::math_utils::atan(mTrackletsCalib[trkWork.getTrackletIndex(iLayer)].getDy() / Geometry::cdrHght()) * TMath::RadToDeg(); + int trkltId = trkWork.getTrackletIndex(iLayer); + // tracklet angle, corrected for pad tilt + const PadPlane* pad = Geometry::instance()->getPadPlane(mTrackletsRaw[trkltId].getDetector()); + float tilt = tan(TMath::DegToRad() * pad->getTiltingAngle()); // tilt is signed! and returned in degrees + float tiltCorrUp = tilt * trkWork.getTgl() * Geometry::cdrHght(); + float padLength = pad->getRowSize(mTrackletsRaw[trkltId].getPadRow()); + if (!((trkWork.getSigmaZ2() < (padLength * padLength / 12.f)) && (std::fabs(mTrackletsCalib[trkltId].getZ() - trkWork.getZ()) < padLength))) { + tiltCorrUp = 0.f; + } + // use uncalibrated dy because online calibration does not work otherwise + float trkltDy = mTrackletsRaw[trkltId].getUncalibratedDy(30.f / o2::trd::constants::VDRIFTDEFAULT) + tiltCorrUp; + float trkltAngle = o2::math_utils::atan(trkltDy / Geometry::cdrHght()) * TMath::RadToDeg(); float angleDeviation = trkltAngle - trkAngle; - if (mAngResHistos.addEntry(angleDeviation, trkAngle, mTrackletsRaw[trkWork.getTrackletIndex(iLayer)].getDetector())) { + if (mAngResHistos.addEntry(angleDeviation, trkAngle, mTrackletsRaw[trkltId].getDetector())) { // track impact angle out of histogram range continue; } diff --git a/Detectors/TRD/workflow/src/TRDTrackletTransformerSpec.cxx b/Detectors/TRD/workflow/src/TRDTrackletTransformerSpec.cxx index 77668ab96a9fd..65826fa76609d 100644 --- a/Detectors/TRD/workflow/src/TRDTrackletTransformerSpec.cxx +++ b/Detectors/TRD/workflow/src/TRDTrackletTransformerSpec.cxx @@ -141,7 +141,13 @@ void TRDTrackletTransformerSpec::finaliseCCDB(ConcreteDataMatcher& matcher, void } if (matcher == ConcreteDataMatcher("TRD", "CALVDRIFTEXB", 0)) { LOG(info) << "CalVdriftExB object has been updated"; - mTransformer.setCalVdriftExB((const o2::trd::CalVdriftExB*)obj); + //for (int iDet = 0; iDet < 540; iDet++) LOGP(info, "vdexb values: {} {} {}", ((o2::trd::CalVdriftExB*)obj)->getVdriftDefaultAvg(iDet), ((o2::trd::CalVdriftExB*)obj)->getAverageVdrift(), ((o2::trd::CalVdriftExB*)obj)->isGoodVdrift(iDet)); + for (int iDet = 0; iDet < constants::MAXCHAMBER; iDet++) { + // set to average value if the calibration is not correct + mTransformer.setVdrift(iDet, ((const o2::trd::CalVdriftExB*)obj)->getVdrift(iDet, true)); + mTransformer.setExB(iDet, ((const o2::trd::CalVdriftExB*)obj)->getExB(iDet, true)); + } + //mTransformer.setCalVdriftExB((const o2::trd::CalVdriftExB*)obj); return; } } diff --git a/GPU/GPUTracking/Definitions/GPUSettingsList.h b/GPU/GPUTracking/Definitions/GPUSettingsList.h index d70fac115eab7..7c2829b95ce19 100644 --- a/GPU/GPUTracking/Definitions/GPUSettingsList.h +++ b/GPU/GPUTracking/Definitions/GPUSettingsList.h @@ -184,6 +184,7 @@ AddOptionRTC(extraRoadZ, float, 10.f, "", 0, "Addition to search road around tra AddOptionRTC(trkltResRPhiIdeal, float, 1.f, "", 0, "Optimal tracklet rphi resolution in cm (in case phi of track = lorentz angle)") AddOptionRTC(maxChi2Red, float, 99.f, "", 0, "maximum chi2 per attached tracklet for TRD tracks TODO: currently effectively disabled, requires tuning") AddOptionRTC(applyDeflectionCut, uint8_t, 0, "", 0, "Set to 1 to enable tracklet selection based on deflection") +AddOptionRTC(addDeflectionInChi2, uint8_t, 0, "", 0, "Set to 1 to add the deflection in the chi2 calculation for matching") AddOptionRTC(stopTrkAfterNMissLy, uint8_t, 6, "", 0, "Abandon track following after N layers without a TRD match") AddOptionRTC(nTrackletsMin, uint8_t, 3, "", 0, "Tracks with less attached tracklets are discarded after the tracking") AddOptionRTC(matCorrType, uint8_t, 2, "", 0, "Material correction to use: 0 - none, 1 - TGeo, 2 - matLUT") diff --git a/GPU/GPUTracking/TRDTracking/GPUTRDTracker.cxx b/GPU/GPUTracking/TRDTracking/GPUTRDTracker.cxx index 2f754d2416bc1..d93b3ff0275ef 100644 --- a/GPU/GPUTracking/TRDTracking/GPUTRDTracker.cxx +++ b/GPU/GPUTracking/TRDTracking/GPUTRDTracker.cxx @@ -92,7 +92,7 @@ void* GPUTRDTracker_t::SetPointersTracks(void* base) } template -GPUTRDTracker_t::GPUTRDTracker_t() : mR(nullptr), mIsInitialized(false), mGenerateSpacePoints(false), mProcessPerTimeFrame(false), mNAngleHistogramBins(25), mAngleHistogramRange(50), mMemoryPermanent(-1), mMemoryTracklets(-1), mMemoryTracks(-1), mNMaxCollisions(0), mNMaxTracks(0), mNMaxSpacePoints(0), mTracks(nullptr), mTrackAttribs(nullptr), mNCandidates(1), mNTracks(0), mNEvents(0), mMaxBackendThreads(100), mTrackletIndexArray(nullptr), mHypothesis(nullptr), mCandidates(nullptr), mSpacePoints(nullptr), mGeo(nullptr), mRPhiA2(0), mRPhiB(0), mRPhiC2(0), mDyA2(0), mDyB(0), mDyC2(0), mAngleToDyA(0), mAngleToDyB(0), mAngleToDyC(0), mDebugOutput(false), mMaxEta(0.84f), mRoadZ(18.f), mZCorrCoefNRC(1.4f), mTPCVdrift(2.58f), mTPCTDriftOffset(0.f), mDebug(new GPUTRDTrackerDebug()) +GPUTRDTracker_t::GPUTRDTracker_t() : mR(nullptr), mIsInitialized(false), mGenerateSpacePoints(false), mProcessPerTimeFrame(false), mNAngleHistogramBins(25), mAngleHistogramRange(50), mMemoryPermanent(-1), mMemoryTracklets(-1), mMemoryTracks(-1), mNMaxCollisions(0), mNMaxTracks(0), mNMaxSpacePoints(0), mTracks(nullptr), mTrackAttribs(nullptr), mNCandidates(1), mNTracks(0), mNEvents(0), mMaxBackendThreads(100), mTrackletIndexArray(nullptr), mHypothesis(nullptr), mCandidates(nullptr), mSpacePoints(nullptr), mGeo(nullptr), mRPhiA2(0), mRPhiB(0), mRPhiC2(0), mDyA2(0), mDyB(0), mDyC2(0), mAngleToDyA(0), mAngleToDyB(0), mAngleToDyC(0), mCorrYDy(0), mDebugOutput(false), mMaxEta(0.84f), mRoadZ(18.f), mZCorrCoefNRC(1.4f), mTPCVdrift(2.58f), mTPCTDriftOffset(0.f), mDebug(new GPUTRDTrackerDebug()) { //-------------------------------------------------------------------- // Default constructor @@ -139,7 +139,7 @@ void GPUTRDTracker_t::UpdateGeometry() if (Bz > 0) { GPUInfo("Loading error parameterization for Bz = +2 kG"); mRPhiA2 = resRPhiIdeal2, mRPhiB = -1.43e-2f, mRPhiC2 = 4.55e-2f; - mDyA2 = 1.225e-3f, mDyB = -9.8e-3f, mDyC2 = 3.88e-2f; + mDyA2 = 1.225e-3f, mDyB = -9.8e-3f, mDyC2 = 3.88e-2f; //TODO: update param also for low field mAngleToDyA = -0.1f, mAngleToDyB = 1.89f, mAngleToDyC = -0.4f; } else { GPUInfo("Loading error parameterization for Bz = -2 kG"); @@ -152,12 +152,14 @@ void GPUTRDTracker_t::UpdateGeometry() if (Bz > 0) { GPUInfo("Loading error parameterization for Bz = +5 kG"); mRPhiA2 = resRPhiIdeal2, mRPhiB = 0.125f, mRPhiC2 = 0.0961f; - mDyA2 = 1.681e-3f, mDyB = 0.15f, mDyC2 = 0.1849f; + //mDyA2 = 1.681e-3f, mDyB = 0.15f, mDyC2 = 0.1849f; + mDyA2 = 7.8e-3f, mDyB = 0.11f, mDyC2 = 0.29f; mAngleToDyA = 0.13f, mAngleToDyB = 2.43f, mAngleToDyC = -0.58f; } else { GPUInfo("Loading error parameterization for Bz = -5 kG"); mRPhiA2 = resRPhiIdeal2, mRPhiB = -0.14f, mRPhiC2 = 0.1156f; - mDyA2 = 2.209e-3f, mDyB = -0.15f, mDyC2 = 0.2025f; + //mDyA2 = 2.209e-3f, mDyB = -0.15f, mDyC2 = 0.2025f; + mDyA2 = 9.1e-3f, mDyB = -0.14f, mDyC2 = 0.35f; mAngleToDyA = -0.15f, mAngleToDyB = 2.34f, mAngleToDyC = 0.56f; } } else { @@ -166,6 +168,8 @@ void GPUTRDTracker_t::UpdateGeometry() GPUWarning("No error parameterization available for Bz = %.2f kG. Keeping default value (sigma_y = const. = 1cm)", Bz); mRPhiA2 = 1.f; } + // Covariance matrix with dy used in chi2 + mCorrYDy = 0.15; // TODO: better understanding of this parameter + dependence on phi, run, bz, ... // obtain average radius of TRD chambers float x0[kNLayers] = {300.2f, 312.8f, 325.4f, 338.0f, 350.6f, 363.2f}; // used as default value in case no transformation matrix can be obtained @@ -243,6 +247,7 @@ void GPUTRDTracker_t::PrepareTracking(GPUChainTracking* chainTrack chainTracking->mIOPtrs.trdSpacePoints = mSpacePoints; } mNEvents++; + } template @@ -615,9 +620,11 @@ GPUd() bool GPUTRDTracker_t::FollowProlongation(PROP* prop, TRDTRK prop->getPropagatedYZ(spacePoints[trkltIdx].getX(), projY, projZ); // correction for tilted pads (only applied if deltaZ < lPad && track z err << lPad) float tiltCorr = tilt * (spacePoints[trkltIdx].getZ() - projZ); + float dyTiltCorr = tilt * trkWork->getTgl() * mGeo->GetCdrHght(); float lPad = pad->GetRowSize(tracklets[trkltIdx].GetZbin()); if (!((CAMath::Abs(spacePoints[trkltIdx].getZ() - projZ) < lPad) && (trkWork->getSigmaZ2() < (lPad * lPad / 12.f)))) { tiltCorr = 0.f; // will be zero also for TPC tracks which are shifted in z + dyTiltCorr = 0.f; } // correction for mean z position of tracklet (is not the center of the pad if track eta != 0) float zPosCorr = spacePoints[trkltIdx].getZ() + mZCorrCoefNRC * trkWork->getTgl(); @@ -626,13 +633,23 @@ GPUd() bool GPUTRDTracker_t::FollowProlongation(PROP* prop, TRDTRK float deltaY = yPosCorr - projY; float deltaZ = zPosCorr - projZ; float trkltPosTmpYZ[2] = {yPosCorr, zPosCorr}; - float trkltCovTmp[3] = {0.f}; + float trkltCovTmpWithDy[6] = {0.f}; if ((CAMath::Abs(deltaY) < roadY) && (CAMath::Abs(deltaZ) < roadZ)) { // TODO: check if this is still necessary after the cut before propagation of track // tracklet is in windwow: get predicted chi2 for update and store tracklet index if best guess - RecalcTrkltCov(tilt, trkWork->getSnp(), pad->GetRowSize(tracklets[trkltIdx].GetZbin()), trkltCovTmp); + RecalcTrkltCov(tilt, trkWork->getSnp(), pad->GetRowSize(tracklets[trkltIdx].GetZbin()), trkltCovTmpWithDy, Param().rec.trd.addDeflectionInChi2); + float trkltCovTmp[3] = {trkltCovTmpWithDy[0], trkltCovTmpWithDy[1], trkltCovTmpWithDy[2]}; float chi2 = prop->getPredictedChi2(trkltPosTmpYZ, trkltCovTmp); + if (Param().rec.trd.addDeflectionInChi2 && (trkWork->getSnp() < 1.f - 1e-6f) && (trkWork->getSnp() > -1.f + 1e-6f)) { + // we add the slope in the chi2 calculation + if (InvertCov(trkltCovTmpWithDy)) { + float deltaDy = spacePoints[trkltIdx].getDy() + dyTiltCorr - ConvertAngleToDy(trkWork->getSnp()); + chi2 = deltaY * trkltCovTmpWithDy[0] * deltaY + 2 * deltaY * trkltCovTmpWithDy[1] * deltaZ + 2 * deltaY * trkltCovTmpWithDy[3] * deltaDy + + deltaZ * trkltCovTmpWithDy[2] * deltaZ + 2 * deltaZ * trkltCovTmpWithDy[4] * deltaDy + + deltaDy * trkltCovTmpWithDy[5] * deltaDy; + } + } // TODO cut on angular pull should be made stricter when proper v-drift calibration for the TRD tracklets is implemented - if ((chi2 > Param().rec.trd.maxChi2) || (Param().rec.trd.applyDeflectionCut && CAMath::Abs(GetAngularPull(spacePoints[trkltIdx].getDy(), trkWork->getSnp())) > 4)) { + if ((chi2 > Param().rec.trd.maxChi2) || (Param().rec.trd.applyDeflectionCut && CAMath::Abs(GetAngularPull(spacePoints[trkltIdx].getDy() + dyTiltCorr, trkWork->getSnp())) > 4)) { continue; } Hypothesis hypo(trkWork->getNlayersFindable(), iCandidate, trkltIdx, trkWork->getChi2() + chi2); @@ -717,9 +734,10 @@ GPUd() bool GPUTRDTracker_t::FollowProlongation(PROP* prop, TRDTRK tiltCorrUp = 0.f; } float trkltPosUp[2] = {spacePoints[mHypothesis[iUpdate + hypothesisIdxOffset].mTrackletId].getY() - tiltCorrUp, zPosCorrUp}; - float trkltCovUp[3] = {0.f}; - RecalcTrkltCov(tilt, trkWork->getSnp(), pad->GetRowSize(tracklets[mHypothesis[iUpdate + hypothesisIdxOffset].mTrackletId].GetZbin()), trkltCovUp); - + float trkltCovUpWithDy[6] = {0.f}; + RecalcTrkltCov(tilt, trkWork->getSnp(), pad->GetRowSize(tracklets[mHypothesis[iUpdate + hypothesisIdxOffset].mTrackletId].GetZbin()), trkltCovUpWithDy, false); + float trkltCovUp[3] = {trkltCovUpWithDy[0], trkltCovUpWithDy[1], trkltCovUpWithDy[2]}; + #ifdef ENABLE_GPUTRDDEBUG prop->setTrack(&trackNoUp); prop->rotate(GetAlphaOfSector(trkltSec)); @@ -959,7 +977,7 @@ GPUd() float GPUTRDTracker_t::GetAlphaOfSector(const int32_t sec) } template -GPUd() void GPUTRDTracker_t::RecalcTrkltCov(const float tilt, const float snp, const float rowSize, float (&cov)[3]) +GPUd() void GPUTRDTracker_t::RecalcTrkltCov(const float tilt, const float snp, const float rowSize, float (&cov)[6], bool withDy) { //-------------------------------------------------------------------- // recalculate tracklet covariance taking track phi angle into account @@ -972,11 +990,71 @@ GPUd() void GPUTRDTracker_t::RecalcTrkltCov(const float tilt, cons cov[0] = c2 * (sy2 + t2 * sz2); cov[1] = c2 * tilt * (sz2 - sy2); cov[2] = c2 * (t2 * sy2 + sz2); + + if (withDy) { + float sdy2 = GetAngularResolution(snp); + cov[3] = mCorrYDy * CAMath::Sqrt(sdy2 * c2 * sy2); + cov[4] = -tilt * mCorrYDy * CAMath::Sqrt(sdy2 * c2 * sy2); + cov[5] = sdy2; + } +} + +template +GPUd() bool GPUTRDTracker_t::InvertCov(float (&cov)[6]) { + // invert a 3*3 symmetric matrix. Adapted from https://root.cern.ch/doc/master/TMatrixTSymCramerInv_8cxx_source.html + + float c00 = cov[2] * cov[5] - cov[4] * cov[4]; + float c01 = cov[4] * cov[3] - cov[1] * cov[5]; + float c02 = cov[1] * cov[4] - cov[2] * cov[3]; + float c11 = cov[5] * cov[0] - cov[3] * cov[3]; + float c12 = cov[3] * cov[1] - cov[4] * cov[0]; + float c22 = cov[0] * cov[2] - cov[1] * cov[1]; + + float t0 = CAMath::Abs(cov[0]); + float t1 = CAMath::Abs(cov[1]); + float t2 = CAMath::Abs(cov[3]); + + float det; + float tmp; + + if (t0 >= t1) { + if (t2 >= t0) { + tmp = cov[3]; + det = c12*c01-c11*c02; + } else { + tmp = cov[0]; + det = c11*c22-c12*c12; + } + } else if (t2 >= t1) { + tmp = cov[3]; + det = c12*c01-c11*c02; + } else { + tmp = cov[1]; + det = c02*c12-c01*c22; + } + + if ( det == 0 || tmp == 0) { + return false; + } + + float s = tmp/det; + + cov[0] = s*c00; + cov[1] = s*c01; + cov[3] = s*c02; + cov[2] = s*c11; + cov[4] = s*c12; + cov[5] = s*c22; + + return true; } template GPUd() float GPUTRDTracker_t::GetAngularPull(float dYtracklet, float snp) const { + if (snp > 1.f - 1e-6f || snp < -1.f + 1e-6f) { + return 999.f; + } float dYtrack = ConvertAngleToDy(snp); float dYresolution = GetAngularResolution(snp); if (dYresolution < 1e-6f) { diff --git a/GPU/GPUTracking/TRDTracking/GPUTRDTracker.h b/GPU/GPUTracking/TRDTracking/GPUTRDTracker.h index f8fa0342ee62d..75f60aedbfdf5 100644 --- a/GPU/GPUTracking/TRDTracking/GPUTRDTracker.h +++ b/GPU/GPUTracking/TRDTracking/GPUTRDTracker.h @@ -26,6 +26,7 @@ #include "GPULogging.h" #include "GPUTRDInterfaces.h" + #ifndef GPUCA_GPUCODE_DEVICE #include #endif @@ -116,9 +117,10 @@ class GPUTRDTracker_t : public GPUProcessor GPUd() float GetAlphaOfSector(const int32_t sec) const; GPUd() float GetRPhiRes(float snp) const { return (mRPhiA2 + mRPhiC2 * (snp - mRPhiB) * (snp - mRPhiB)); } // parametrization obtained from track-tracklet residuals: GPUd() float GetAngularResolution(float snp) const { return mDyA2 + mDyC2 * (snp - mDyB) * (snp - mDyB); } // a^2 + c^2 * (snp - b)^2 - GPUd() float ConvertAngleToDy(float snp) const { return mAngleToDyA + mAngleToDyB * snp + mAngleToDyC * snp * snp; } // a + b*snp + c*snp^2 is more accurate than sin(phi) = (dy / xDrift) / sqrt(1+(dy/xDrift)^2) + GPUd() float ConvertAngleToDy(float snp) const { return 3.f * snp / sqrt(1 - snp * snp); } // when calibrated, sin(phi) = (dy / xDrift) / sqrt(1+(dy/xDrift)^2) works well GPUd() float GetAngularPull(float dYtracklet, float snp) const; - GPUd() void RecalcTrkltCov(const float tilt, const float snp, const float rowSize, float (&cov)[3]); + GPUd() void RecalcTrkltCov(const float tilt, const float snp, const float rowSize, float (&cov)[6], bool withDy = false); + GPUd() bool InvertCov(float (&cov)[6]); GPUd() void FindChambersInRoad(const TRDTRK* t, const float roadY, const float roadZ, const int32_t iLayer, int32_t* det, const float zMax, const float alpha, const float zShiftTrk) const; GPUd() bool IsGeoFindable(const TRDTRK* t, const int32_t layer, const float alpha, const float zShiftTrk) const; GPUd() void InsertHypothesis(Hypothesis hypo, int32_t& nCurrHypothesis, int32_t idxOffset); @@ -184,6 +186,7 @@ class GPUTRDTracker_t : public GPUProcessor float mAngleToDyA; // parameterization for conversion track angle -> tracklet deflection float mAngleToDyB; // parameterization for conversion track angle -> tracklet deflection float mAngleToDyC; // parameterization for conversion track angle -> tracklet deflection + float mCorrYDy; // correlation coefficient between tracklet position and tracklet deflection /// ---- end error parametrization ---- bool mDebugOutput; // store debug output static constexpr const float sRadialOffset = -0.1f; // due to (possible) mis-calibration of t0 -> will become obsolete when tracklet conversion is done outside of the tracker