diff --git a/DataFormats/Detectors/TRD/include/DataFormatsTRD/CalGain.h b/DataFormats/Detectors/TRD/include/DataFormatsTRD/CalGain.h index 408198665df6a..039151e6890cf 100644 --- a/DataFormats/Detectors/TRD/include/DataFormatsTRD/CalGain.h +++ b/DataFormats/Detectors/TRD/include/DataFormatsTRD/CalGain.h @@ -33,22 +33,25 @@ class CalGain void setMPVdEdx(int iDet, float mpv) { mMPVdEdx[iDet] = mpv; } - float getMPVdEdx(int iDet, bool defaultAvg = false) const { + 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(); + if (!defaultAvg || isGoodGain(iDet)) + return mMPVdEdx[iDet]; + else + return getAverageGain(); } - - - float getAverageGain() const { + + 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 ++; + ngood++; averageGain += mMPVdEdx[iDet]; } } @@ -59,10 +62,13 @@ class CalGain averageGain /= ngood; return averageGain; } - - bool isGoodGain(int iDet) const { - if (TMath::Abs(mMPVdEdx[iDet] - constants::MPVDEDXDEFAULT) > 1e-6) return true; - else return false; + + bool isGoodGain(int iDet) const + { + if (TMath::Abs(mMPVdEdx[iDet] - constants::MPVDEDXDEFAULT) > 1e-6) + return true; + else + return false; } private: diff --git a/DataFormats/Detectors/TRD/include/DataFormatsTRD/CalVdriftExB.h b/DataFormats/Detectors/TRD/include/DataFormatsTRD/CalVdriftExB.h index aa47a75a65e7a..58eba0e959a31 100644 --- a/DataFormats/Detectors/TRD/include/DataFormatsTRD/CalVdriftExB.h +++ b/DataFormats/Detectors/TRD/include/DataFormatsTRD/CalVdriftExB.h @@ -33,26 +33,33 @@ 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 { + + 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(); + 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 getExB(int iDet, bool defaultAvg = false) const + { + if (!defaultAvg || (isGoodExB(iDet) && isGoodVdrift(iDet))) + return mExB[iDet]; + else + return getAverageExB(); } - - float getAverageVdrift() const { + + 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 ++; + ngood++; averageVdrift += mVdrift[iDet]; } } @@ -64,14 +71,15 @@ class CalVdriftExB return averageVdrift; } - float getAverageExB() const { + 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 ++; + ngood++; averageExB += mExB[iDet]; } } @@ -83,7 +91,8 @@ class CalVdriftExB return averageExB; } - bool isGoodExB(int iDet) const { + 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 @@ -91,10 +100,12 @@ class CalVdriftExB TMath::Abs(mExB[iDet] - constants::EXBMIN) > 0.01 && TMath::Abs(mExB[iDet] - constants::EXBMAX) > 0.01) return true; - else return false; + else + return false; } - bool isGoodVdrift(int iDet) const { + 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 @@ -102,7 +113,8 @@ class CalVdriftExB TMath::Abs(mVdrift[iDet] - constants::VDRIFTMIN) > 0.1 && TMath::Abs(mVdrift[iDet] - constants::VDRIFTMAX) > 0.1) return true; - else return false; + else + return false; } private: diff --git a/DataFormats/Detectors/TRD/include/DataFormatsTRD/Constants.h b/DataFormats/Detectors/TRD/include/DataFormatsTRD/Constants.h index b2c827710356d..9a4da1024e251 100644 --- a/DataFormats/Detectors/TRD/include/DataFormatsTRD/Constants.h +++ b/DataFormats/Detectors/TRD/include/DataFormatsTRD/Constants.h @@ -75,11 +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 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 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 c4367ef038cbc..2b85de743e955 100644 --- a/Detectors/Align/Workflow/src/BarrelAlignmentSpec.cxx +++ b/Detectors/Align/Workflow/src/BarrelAlignmentSpec.cxx @@ -312,7 +312,7 @@ void BarrelAlignmentSpec::finaliseCCDB(o2::framework::ConcreteDataMatcher& match 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); + // 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 1defe41ec9f58..a9edad2ea7733 100644 --- a/Detectors/TRD/base/include/TRDBase/TrackletTransformer.h +++ b/Detectors/TRD/base/include/TRDBase/TrackletTransformer.h @@ -30,13 +30,11 @@ class TrackletTransformer void init(); - void setVdrift(int iDet, float vd) {mVdrift[iDet] = vd;} - void setExB(int iDet, float exb) {mExB[iDet] = exb;} + 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; diff --git a/Detectors/TRD/base/src/TrackletTransformer.cxx b/Detectors/TRD/base/src/TrackletTransformer.cxx index 370426d7c7201..08b60a92e3cc7 100644 --- a/Detectors/TRD/base/src/TrackletTransformer.cxx +++ b/Detectors/TRD/base/src/TrackletTransformer.cxx @@ -40,8 +40,8 @@ float TrackletTransformer::calculateDy(int detector, int slope, const PadPlane* { double padWidth = padPlane->getWidthIPad(); - //float vDrift = mCalVdriftExB->getVdrift(detector, true); - //float exb = mCalVdriftExB->getExB(detector, true); + // float vDrift = mCalVdriftExB->getVdrift(detector, true); + // float exb = mCalVdriftExB->getExB(detector, true); float vDrift = mVdrift[detector]; float exb = mExB[detector]; @@ -57,7 +57,7 @@ float TrackletTransformer::calculateDy(int detector, int slope, const PadPlane* // assuming angle in Bailhache, fig. 4.17 would be positive in our calibration code double calibratedDy = rawDy - lorentzCorrection; - + return calibratedDy; } @@ -99,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); diff --git a/Detectors/TRD/calibration/include/TRDCalibration/CalibrationParams.h b/Detectors/TRD/calibration/include/TRDCalibration/CalibrationParams.h index 31f7f56850419..cadd21af2a55d 100644 --- a/Detectors/TRD/calibration/include/TRDCalibration/CalibrationParams.h +++ b/Detectors/TRD/calibration/include/TRDCalibration/CalibrationParams.h @@ -25,9 +25,9 @@ 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 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 = 200; ///< minimum number of entries per chamber to fit single time slot + 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 diff --git a/Detectors/TRD/calibration/include/TRDCalibration/CalibratorVdExB.h b/Detectors/TRD/calibration/include/TRDCalibration/CalibratorVdExB.h index 50cc3930baa96..16223f3e78112 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 c9eb5ec0d66a6..23f6a1dbb008d 100644 --- a/Detectors/TRD/calibration/macros/manualCalibFit.C +++ b/Detectors/TRD/calibration/macros/manualCalibFit.C @@ -60,10 +60,10 @@ void manualCalibFit(int runNumber = 563335, bool usePreCorrFromCCDB = false) mNEntriesPerBinSum.fill(0); tree->SetBranchAddress("mHistogramEntries[13500]", &mHistogramEntries); tree->SetBranchAddress("mNEntriesPerBin[13500]", &mNEntriesPerBin); - - // use precorr values from ccdb + + // 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(); @@ -71,13 +71,12 @@ void manualCalibFit(int runNumber = 563335, bool usePreCorrFromCCDB = false) 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 @@ -88,9 +87,11 @@ void manualCalibFit(int runNumber = 563335, bool usePreCorrFromCCDB = false) 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) << " " << calObject->getExB(iDet) << " "; + if (calObject->isGoodExB(iDet)) + counter++; + if (iDet % 6 == 5) + std::cout << std::endl; mFitFunctor.vdPreCorr[iDet] = calObject->getVdrift(iDet, true); mFitFunctor.laPreCorr[iDet] = calObject->getExB(iDet, true); } @@ -128,7 +129,8 @@ void manualCalibFit(int runNumber = 563335, bool usePreCorrFromCCDB = false) 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)<GetBinEntries(iBin + 1) << " " << mFitFunctor.profiles[iDet]->GetBinEffectiveEntries(iBin + 1) << " " << mFitFunctor.profiles[iDet]->GetBinError(iBin + 1) << std::endl; } } printf("Det %d: nEntries=%d \n", iDet, nEntriesDetTotal[iDet]); @@ -140,13 +142,14 @@ void manualCalibFit(int runNumber = 563335, bool usePreCorrFromCCDB = false) printf("-------- Started fits\n"); std::array 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; + if (nEntriesDetTotal[iDet] < 75) + continue; mFitFunctor.currDet = iDet; ROOT::Fit::Fitter fitter; double paramsStart[2]; @@ -168,19 +171,20 @@ void manualCalibFit(int runNumber = 563335, bool usePreCorrFromCCDB = false) auto fitResult = fitter.Result(); laFitResults[iDet] = fitResult.Parameter(0); vdFitResults[iDet] = fitResult.Parameter(1); - 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]); + 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()); + 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()<GetEntries() << std::endl; + ; + std::cout << "vdrift mean: " << hVd->GetMean() << " sigma: " << hVd->GetStdDev() << std::endl; + std::cout << "lorentz angle mean: " << hLa->GetMean() << " sigma: " << hLa->GetStdDev() << std::endl; //---------------------------------------------------- // Write @@ -191,5 +195,4 @@ void manualCalibFit(int runNumber = 563335, bool usePreCorrFromCCDB = false) 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 68f2998dd26a2..7aae6d14d5023 100644 --- a/Detectors/TRD/calibration/src/CalibratorVdExB.cxx +++ b/Detectors/TRD/calibration/src/CalibratorVdExB.cxx @@ -24,7 +24,7 @@ #include "CCDB/BasicCCDBManager.h" #include "CommonUtils/NameConf.h" #include "CommonUtils/MemFileHelper.h" -//#include "DetectorsBase/Propagator.h" +// #include "DetectorsBase/Propagator.h" #include #include @@ -109,14 +109,14 @@ void CalibratorVdExB::initProcessing() // 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(); + // 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(); diff --git a/Detectors/TRD/calibration/src/TrackBasedCalib.cxx b/Detectors/TRD/calibration/src/TrackBasedCalib.cxx index e794e14873d1f..0d87bdb2b3014 100644 --- a/Detectors/TRD/calibration/src/TrackBasedCalib.cxx +++ b/Detectors/TRD/calibration/src/TrackBasedCalib.cxx @@ -187,8 +187,8 @@ 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 - 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 + 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; } } diff --git a/Detectors/TRD/workflow/src/TRDTrackletTransformerSpec.cxx b/Detectors/TRD/workflow/src/TRDTrackletTransformerSpec.cxx index 65826fa76609d..038d1c87c1d49 100644 --- a/Detectors/TRD/workflow/src/TRDTrackletTransformerSpec.cxx +++ b/Detectors/TRD/workflow/src/TRDTrackletTransformerSpec.cxx @@ -141,13 +141,13 @@ void TRDTrackletTransformerSpec::finaliseCCDB(ConcreteDataMatcher& matcher, void } if (matcher == ConcreteDataMatcher("TRD", "CALVDRIFTEXB", 0)) { LOG(info) << "CalVdriftExB object has been updated"; - //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 < 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); + // mTransformer.setCalVdriftExB((const o2::trd::CalVdriftExB*)obj); return; } } diff --git a/GPU/GPUTracking/TRDTracking/GPUTRDTracker.cxx b/GPU/GPUTracking/TRDTracking/GPUTRDTracker.cxx index d93b3ff0275ef..51e37471a4881 100644 --- a/GPU/GPUTracking/TRDTracking/GPUTRDTracker.cxx +++ b/GPU/GPUTracking/TRDTracking/GPUTRDTracker.cxx @@ -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; //TODO: update param also for low field + 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,13 +152,13 @@ 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; } @@ -247,7 +247,6 @@ void GPUTRDTracker_t::PrepareTracking(GPUChainTracking* chainTrack chainTracking->mIOPtrs.trdSpacePoints = mSpacePoints; } mNEvents++; - } template @@ -643,9 +642,7 @@ GPUd() bool GPUTRDTracker_t::FollowProlongation(PROP* prop, TRDTRK // 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; + 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 @@ -737,7 +734,7 @@ GPUd() bool GPUTRDTracker_t::FollowProlongation(PROP* prop, TRDTRK 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)); @@ -990,7 +987,7 @@ 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); @@ -1000,52 +997,53 @@ GPUd() void GPUTRDTracker_t::RecalcTrkltCov(const float tilt, cons } template -GPUd() bool GPUTRDTracker_t::InvertCov(float (&cov)[6]) { +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 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; + det = c12 * c01 - c11 * c02; } else { tmp = cov[0]; - det = c11*c22-c12*c12; + det = c11 * c22 - c12 * c12; } } else if (t2 >= t1) { tmp = cov[3]; - det = c12*c01-c11*c02; + det = c12 * c01 - c11 * c02; } else { tmp = cov[1]; - det = c02*c12-c01*c22; + det = c02 * c12 - c01 * c22; } - - if ( det == 0 || tmp == 0) { + + 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; - + + 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; } diff --git a/GPU/GPUTracking/TRDTracking/GPUTRDTracker.h b/GPU/GPUTracking/TRDTracking/GPUTRDTracker.h index 75f60aedbfdf5..5d80fb7c01518 100644 --- a/GPU/GPUTracking/TRDTracking/GPUTRDTracker.h +++ b/GPU/GPUTracking/TRDTracking/GPUTRDTracker.h @@ -26,7 +26,6 @@ #include "GPULogging.h" #include "GPUTRDInterfaces.h" - #ifndef GPUCA_GPUCODE_DEVICE #include #endif @@ -117,7 +116,7 @@ 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 3.f * snp / sqrt(1 - snp * snp); } // when calibrated, sin(phi) = (dy / xDrift) / sqrt(1+(dy/xDrift)^2) works well + 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)[6], bool withDy = false); GPUd() bool InvertCov(float (&cov)[6]);