Skip to content

Commit a449f99

Browse files
glegrasglegras
committed
updates vd exb calib
1 parent aab079f commit a449f99

File tree

11 files changed

+252
-35
lines changed

11 files changed

+252
-35
lines changed

‎DataFormats/Detectors/TRD/include/DataFormatsTRD/CalGain.h‎

Lines changed: 32 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -33,12 +33,42 @@ class CalGain
3333

3434
void setMPVdEdx(int iDet, float mpv) { mMPVdEdx[iDet] = mpv; }
3535

36-
float getMPVdEdx(int iDet) const { return mMPVdEdx[iDet]; }
36+
float getMPVdEdx(int iDet, bool defaultAvg = false) const {
37+
// if defaultAvg = false, we take the value stored whatever it is
38+
// if defaultAvg = true and we have default value or bad value stored, we take the average on all chambers instead
39+
if (!defaultAvg || isGoodGain(iDet)) return mMPVdEdx[iDet];
40+
else return getAverageGain();
41+
}
42+
43+
44+
float getAverageGain() const {
45+
float averageGain = 0.;
46+
int ngood = 0;
47+
48+
for (int iDet = 0; iDet < constants::MAXCHAMBER; iDet++) {
49+
if (isGoodGain(iDet)) {
50+
// The chamber has correct calibration
51+
ngood ++;
52+
averageGain += mMPVdEdx[iDet];
53+
}
54+
}
55+
if (ngood == 0) {
56+
// we should make sure it never happens
57+
return constants::MPVDEDXDEFAULT;
58+
}
59+
averageGain /= ngood;
60+
return averageGain;
61+
}
62+
63+
bool isGoodGain(int iDet) const {
64+
if (TMath::Abs(mMPVdEdx[iDet] - constants::MPVDEDXDEFAULT) > 1e-6) return true;
65+
else return false;
66+
}
3767

3868
private:
3969
std::array<float, constants::MAXCHAMBER> mMPVdEdx{}; ///< Most probable value of dEdx distribution per TRD chamber
4070

41-
ClassDefNV(CalGain, 1);
71+
ClassDefNV(CalGain, 2);
4272
};
4373

4474
} // namespace trd

‎DataFormats/Detectors/TRD/include/DataFormatsTRD/CalVdriftExB.h‎

Lines changed: 71 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -33,15 +33,83 @@ class CalVdriftExB
3333

3434
void setVdrift(int iDet, float vd) { mVdrift[iDet] = vd; }
3535
void setExB(int iDet, float exb) { mExB[iDet] = exb; }
36+
37+
float getVdrift(int iDet, bool defaultAvg = false) const {
38+
// if defaultAvg = false, we take the value stored whatever it is
39+
// if defaultAvg = true and we have default value or bad value stored, we take the average on all chambers instead
40+
if (!defaultAvg || (isGoodExB(iDet) && isGoodVdrift(iDet))) return mVdrift[iDet];
41+
else return getAverageVdrift();
42+
}
43+
float getExB(int iDet, bool defaultAvg = false) const {
44+
if (!defaultAvg || (isGoodExB(iDet) && isGoodVdrift(iDet))) return mExB[iDet];
45+
else return getAverageExB();
46+
}
47+
48+
float getAverageVdrift() const {
49+
float averageVdrift = 0.;
50+
int ngood = 0;
51+
52+
for (int iDet = 0; iDet < constants::MAXCHAMBER; iDet++) {
53+
if (isGoodExB(iDet) && isGoodVdrift(iDet)) {
54+
// Both values need to be correct to declare a chamber as well calibrated
55+
ngood ++;
56+
averageVdrift += mVdrift[iDet];
57+
}
58+
}
59+
if (ngood == 0) {
60+
// we should make sure it never happens
61+
return constants::VDRIFTDEFAULT;
62+
}
63+
averageVdrift /= ngood;
64+
return averageVdrift;
65+
}
3666

37-
float getVdrift(int iDet) const { return mVdrift[iDet]; }
38-
float getExB(int iDet) const { return mExB[iDet]; }
67+
float getAverageExB() const {
68+
float averageExB = 0.;
69+
int ngood = 0;
70+
71+
for (int iDet = 0; iDet < constants::MAXCHAMBER; iDet++) {
72+
if (isGoodExB(iDet) && isGoodVdrift(iDet)) {
73+
// Both values need to be correct to declare a chamber as well calibrated
74+
ngood ++;
75+
averageExB += mExB[iDet];
76+
}
77+
}
78+
if (ngood == 0) {
79+
// we should make sure it never happens
80+
return constants::EXBDEFAULT;
81+
}
82+
averageExB /= ngood;
83+
return averageExB;
84+
}
85+
86+
bool isGoodExB(int iDet) const {
87+
// check if value is well calibrated or not
88+
// default calibration if not enough entries
89+
// close to boundaries indicate a failed fit
90+
if (TMath::Abs(mExB[iDet] - constants::EXBDEFAULT) > 1e-6 &&
91+
TMath::Abs(mExB[iDet] - constants::EXBMIN) > 0.01 &&
92+
TMath::Abs(mExB[iDet] - constants::EXBMAX) > 0.01)
93+
return true;
94+
else return false;
95+
}
96+
97+
bool isGoodVdrift(int iDet) const {
98+
// check if value is well calibrated or not
99+
// default calibration if not enough entries
100+
// close to boundaries indicate a failed fit
101+
if (TMath::Abs(mVdrift[iDet] - constants::VDRIFTDEFAULT) > 1e-6 &&
102+
TMath::Abs(mVdrift[iDet] - constants::VDRIFTMIN) > 0.1 &&
103+
TMath::Abs(mVdrift[iDet] - constants::VDRIFTMAX) > 0.1)
104+
return true;
105+
else return false;
106+
}
39107

40108
private:
41109
std::array<float, constants::MAXCHAMBER> mVdrift{}; ///< calibrated drift velocity per TRD chamber
42110
std::array<float, constants::MAXCHAMBER> mExB{}; ///< calibrated Lorentz angle per TRD chamber
43111

44-
ClassDefNV(CalVdriftExB, 1);
112+
ClassDefNV(CalVdriftExB, 2);
45113
};
46114

47115
} // namespace trd

‎DataFormats/Detectors/TRD/include/DataFormatsTRD/Constants.h‎

Lines changed: 4 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -75,7 +75,11 @@ constexpr int TIMEBINS = 30; ///< the number of time bins
7575
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
7676
constexpr int NBINSANGLEDIFF = 25; ///< the number of bins for the track angle used for the vDrift and ExB calibration based on the tracking
7777
constexpr double VDRIFTDEFAULT = 1.546; ///< default value for vDrift
78+
constexpr double VDRIFTMIN = 0.01; ///< min value for vDrift
79+
constexpr double VDRIFTMAX = 3.0; ///< max value for vDrift
7880
constexpr double EXBDEFAULT = 0.0; ///< default value for LorentzAngle
81+
constexpr double EXBMIN = -0.7; ///< min value for LorentzAngle
82+
constexpr double EXBMAX = 0.7; ///< max value for LorentzAngle
7983
constexpr int NBINSGAINCALIB = 320; ///< number of bins in the charge (Q0+Q1+Q2) histogram for gain calibration
8084
constexpr float MPVDEDXDEFAULT = 42.; ///< default Most Probable Value of TRD dEdx
8185
constexpr float T0DEFAULT = 1.2; ///< default value for t0

‎Detectors/Align/Workflow/src/BarrelAlignmentSpec.cxx‎

Lines changed: 6 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -307,7 +307,12 @@ void BarrelAlignmentSpec::finaliseCCDB(o2::framework::ConcreteDataMatcher& match
307307
}
308308
if (matcher == ConcreteDataMatcher("TRD", "CALVDRIFTEXB", 0)) {
309309
LOG(info) << "CalVdriftExB object has been updated";
310-
mTRDTransformer->setCalVdriftExB((const o2::trd::CalVdriftExB*)obj);
310+
for (int iDet = 0; iDet < o2::trd::constants::MAXCHAMBER; iDet++) {
311+
// set to average value if the calibration is not correct
312+
mTRDTransformer->setVdrift(iDet, ((const o2::trd::CalVdriftExB*)obj)->getVdrift(iDet, true));
313+
mTRDTransformer->setExB(iDet, ((const o2::trd::CalVdriftExB*)obj)->getExB(iDet, true));
314+
}
315+
//mTRDTransformer->setCalVdriftExB((const o2::trd::CalVdriftExB*)obj);
311316
return;
312317
}
313318
if (mTPCVDriftHelper.accountCCDBInputs(matcher, obj)) {

‎Detectors/TRD/base/include/TRDBase/TrackletTransformer.h‎

Lines changed: 6 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -30,10 +30,13 @@ class TrackletTransformer
3030

3131
void init();
3232

33-
void setCalVdriftExB(const CalVdriftExB* cal) { mCalVdriftExB = cal; };
33+
void setVdrift(int iDet, float vd) {mVdrift[iDet] = vd;}
34+
void setExB(int iDet, float exb) {mExB[iDet] = exb;}
3435
void setApplyXOR() { mApplyXOR = true; }
3536
void setApplyShift(bool f) { mApplyShift = f; }
3637
bool isShiftApplied() const { return mApplyShift; }
38+
39+
3740

3841
float calculateZ(int padrow, const PadPlane* padPlane) const;
3942

@@ -54,7 +57,8 @@ class TrackletTransformer
5457

5558
float mXAnode;
5659

57-
const CalVdriftExB* mCalVdriftExB{nullptr};
60+
std::array<float, constants::MAXCHAMBER> mVdrift{};
61+
std::array<float, constants::MAXCHAMBER> mExB{};
5862
};
5963

6064
} // namespace trd

‎Detectors/TRD/base/src/TrackletTransformer.cxx‎

Lines changed: 9 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -40,8 +40,10 @@ float TrackletTransformer::calculateDy(int detector, int slope, const PadPlane*
4040
{
4141
double padWidth = padPlane->getWidthIPad();
4242

43-
float vDrift = mCalVdriftExB->getVdrift(detector);
44-
float exb = mCalVdriftExB->getExB(detector);
43+
//float vDrift = mCalVdriftExB->getVdrift(detector, true);
44+
//float exb = mCalVdriftExB->getExB(detector, true);
45+
float vDrift = mVdrift[detector];
46+
float exb = mExB[detector];
4547

4648
// dy = slope * nTimeBins * padWidth * GRANULARITYTRKLSLOPE;
4749
// 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*
5052
// NOTE: check what drift height is used in calibration code to ensure consistency
5153
// NOTE: check sign convention of Lorentz angle
5254
// NOTE: confirm the direction in which vDrift is measured/determined. Is it in x or in direction of drift?
53-
double lorentzCorrection = TMath::Tan(exb) * mXAnode;
55+
// The Lorentz correction have to be applied both at the point of entrance and at the end of the drift region
56+
double lorentzCorrection = TMath::Tan(exb) * mGeo->cdrHght();
5457

5558
// assuming angle in Bailhache, fig. 4.17 would be positive in our calibration code
5659
double calibratedDy = rawDy - lorentzCorrection;
57-
60+
5861
return calibratedDy;
5962
}
6063

@@ -96,7 +99,7 @@ CalibratedTracklet TrackletTransformer::transformTracklet(Tracklet64 tracklet, b
9699
position = tracklet.getPositionBinSigned();
97100
slope = tracklet.getSlopeBinSigned();
98101
}
99-
102+
100103
// calculate raw local chamber space point
101104
const auto padPlane = mGeo->getPadPlane(detector);
102105

@@ -124,7 +127,7 @@ CalibratedTracklet TrackletTransformer::transformTracklet(Tracklet64 tracklet, b
124127
double TrackletTransformer::getTimebin(int detector, double x) const
125128
{
126129
// calculate timebin from x position within chamber
127-
float vDrift = mCalVdriftExB->getVdrift(detector);
130+
float vDrift = mVdrift[detector];
128131
double t0 = 4.0; // time (in timebins) of start of drift region
129132

130133
double timebin;

‎Detectors/TRD/calibration/include/TRDCalibration/CalibrationParams.h‎

Lines changed: 4 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -25,12 +25,13 @@ namespace trd
2525
/// VDrift and ExB calibration parameters.
2626
struct TRDCalibParams : public o2::conf::ConfigurableParamHelper<TRDCalibParams> {
2727
unsigned int nTrackletsMin = 5; ///< minimum amount of tracklets
28+
unsigned int nTrackletsMinLoose = 4; ///< minimum amount of tracklets if two layers with a large lever arm both have a hit
2829
unsigned int chi2RedMax = 6; ///< maximum reduced chi2 acceptable for track quality
29-
size_t minEntriesChamber = 75; ///< minimum number of entries per chamber to fit single time slot
30-
size_t minEntriesTotal = 40'500; ///< minimum total required for meaningful fits
30+
size_t minEntriesChamber = 200; ///< minimum number of entries per chamber to fit single time slot
31+
size_t minEntriesTotal = 400'000; ///< minimum total required for meaningful fits
3132

3233
// For gain calibration
33-
unsigned int nTrackletsMinGainCalib = 5;
34+
unsigned int nTrackletsMinGainCalib = 3;
3435
size_t minEntriesChamberGainCalib = 500; ///< minimum number of entries per chamber to fit single time slot
3536
size_t minEntriesTotalGainCalib = 1'000'000; ///< minimum total required for meaningful fits
3637
// Cuts for selecting clean pion candidates for gain calibration

‎Detectors/TRD/calibration/macros/manualCalibFit.C‎

Lines changed: 65 additions & 9 deletions
Original file line numberDiff line numberDiff line change
@@ -30,18 +30,19 @@
3030

3131
// O2 header
3232
#include <TRDCalibration/CalibratorVdExB.h>
33+
#include "DetectorsBase/Propagator.h"
3334

3435
#endif
3536

3637
// This root macro reads in 'trdangreshistos.root' and
3738
// performs the calibration fits manually as in CalibratorVdExB.cxx
3839
// This can be used for checking if the calibration fits make sense.
39-
void manualCalibFit()
40+
void manualCalibFit(int runNumber = 563336, bool usePreCorrFromCCDB = true)
4041
{
4142
//----------------------------------------------------
4243
// TTree and File
4344
//----------------------------------------------------
44-
std::unique_ptr<TFile> inFilePtr(TFile::Open("trdangreshistos.root"));
45+
std::unique_ptr<TFile> inFilePtr(TFile::Open("trdcaliboutput.root"));
4546
if (inFilePtr == nullptr) {
4647
printf("Input File could not be read!\n'");
4748
return;
@@ -59,19 +60,49 @@ void manualCalibFit()
5960
mNEntriesPerBinSum.fill(0);
6061
tree->SetBranchAddress("mHistogramEntries[13500]", &mHistogramEntries);
6162
tree->SetBranchAddress("mNEntriesPerBin[13500]", &mNEntriesPerBin);
63+
64+
// use precorr values from ccdb
65+
// necessary when the angular residuals were calculated already using ccdb calibration (e.g. in a local run)
66+
67+
o2::trd::CalVdriftExB* calObject;
68+
if (usePreCorrFromCCDB) {
69+
auto& ccdbmgr = o2::ccdb::BasicCCDBManager::instance();
70+
71+
o2::ccdb::CcdbApi ccdb;
72+
ccdb.init("http://alice-ccdb.cern.ch");
73+
auto runDuration = ccdbmgr.getRunDuration(runNumber);
74+
75+
std::map<std::string, std::string> metadata;
76+
std::map<std::string, std::string> headers;
77+
78+
calObject = ccdb.retrieveFromTFileAny<o2::trd::CalVdriftExB>("TRD/Calib/CalVdriftExB", metadata, runDuration.first + 60000, &headers, "", "", "1689478811721");
79+
}
80+
6281

6382
//----------------------------------------------------
6483
// Configure Fitter
6584
//----------------------------------------------------
6685
o2::trd::FitFunctor mFitFunctor;
6786
std::array<std::unique_ptr<TProfile>, 540> profiles; ///< profile histograms for each TRD chamber
87+
int counter = 0;
6888
for (int iDet = 0; iDet < 540; ++iDet) {
6989
mFitFunctor.profiles[iDet] = std::make_unique<TProfile>(Form("profAngleDiff_%i", iDet), Form("profAngleDiff_%i", iDet), 25, -25.f, 25.f);
90+
if (usePreCorrFromCCDB) {
91+
std::cout<<iDet<<" "<<calObject->getVdrift(iDet)<<" "<<calObject->getExB(iDet)<<" ";
92+
if (calObject->isGoodExB(iDet)) counter++;
93+
if (iDet%6==5)std::cout<<std::endl;
94+
mFitFunctor.vdPreCorr[iDet] = calObject->getVdriftDefaultAvg(iDet);
95+
mFitFunctor.laPreCorr[iDet] = calObject->getExBDefaultAvg(iDet);
96+
}
97+
}
98+
std::cout << counter << " good entries in the CCDB " << std::endl;
99+
mFitFunctor.mAnodePlane = 3.35; // don't really care as long as it's not zero, this parameter could be removed
100+
mFitFunctor.lowerBoundAngleFit = 85 * TMath::DegToRad();
101+
mFitFunctor.upperBoundAngleFit = 105 * TMath::DegToRad();
102+
if (!usePreCorrFromCCDB) {
103+
mFitFunctor.vdPreCorr.fill(1.546);
104+
mFitFunctor.laPreCorr.fill(0.0);
70105
}
71-
mFitFunctor.lowerBoundAngleFit = 80 * TMath::DegToRad();
72-
mFitFunctor.upperBoundAngleFit = 100 * TMath::DegToRad();
73-
mFitFunctor.vdPreCorr.fill(1.546);
74-
mFitFunctor.laPreCorr.fill(0.0);
75106

76107
//----------------------------------------------------
77108
// Loop
@@ -88,15 +119,18 @@ void manualCalibFit()
88119
//----------------------------------------------------
89120
// Fill profiles
90121
//----------------------------------------------------
122+
int nEntriesDetTotal[540] = {};
91123
for (int iDet = 0; iDet < 540; ++iDet) {
92124
for (int iBin = 0; iBin < 25; ++iBin) {
93125
auto angleDiffSum = mHistogramEntriesSum[iDet * 25 + iBin];
94126
auto nEntries = mNEntriesPerBinSum[iDet * 25 + iBin];
127+
nEntriesDetTotal[iDet] += nEntries;
95128
if (nEntries > 0) { // skip entries which have no entries; ?
96129
// add to the respective profile for fitting later on
97130
mFitFunctor.profiles[iDet]->Fill(2 * iBin - 25.f, angleDiffSum / nEntries, nEntries);
98131
}
99132
}
133+
printf("Det %d: nEntries=%d \n", iDet, nEntriesDetTotal[iDet]);
100134
}
101135

102136
//----------------------------------------------------
@@ -105,16 +139,21 @@ void manualCalibFit()
105139
printf("-------- Started fits\n");
106140
std::array<float, 540> laFitResults{};
107141
std::array<float, 540> vdFitResults{};
142+
143+
TH1F* hVd = new TH1F("hVd", "v drift", 150, 0.5, 2.);
144+
TH1F* hLa = new TH1F("hLa", "lorentz angle", 200, -25., 25.);
145+
108146
for (int iDet = 0; iDet < 540; ++iDet) {
147+
if (nEntriesDetTotal[iDet] < 75) continue;
109148
mFitFunctor.currDet = iDet;
110149
ROOT::Fit::Fitter fitter;
111150
double paramsStart[2];
112-
paramsStart[0] = 0. * TMath::DegToRad();
151+
paramsStart[0] = 0.;
113152
paramsStart[1] = 1.;
114153
fitter.SetFCN<o2::trd::FitFunctor>(2, mFitFunctor, paramsStart);
115154
fitter.Config().ParSettings(0).SetLimits(-0.7, 0.7);
116155
fitter.Config().ParSettings(0).SetStepSize(.01);
117-
fitter.Config().ParSettings(1).SetLimits(0., 3.);
156+
fitter.Config().ParSettings(1).SetLimits(0.01, 3.);
118157
fitter.Config().ParSettings(1).SetStepSize(.01);
119158
ROOT::Math::MinimizerOptions opt;
120159
opt.SetMinimizerType("Minuit2");
@@ -127,9 +166,26 @@ void manualCalibFit()
127166
auto fitResult = fitter.Result();
128167
laFitResults[iDet] = fitResult.Parameter(0);
129168
vdFitResults[iDet] = fitResult.Parameter(1);
130-
printf("Det %d: la=%f\tvd=%f\n", iDet, laFitResults[iDet] * TMath::RadToDeg(), vdFitResults[iDet]);
169+
printf("Det %d: la=%.3f +/- %.3f (precorr=:%.3f) \tvd=%.3f +/- %.3f (precorr=:%.3f) \t100*minValue=%f\n", iDet, laFitResults[iDet] * TMath::RadToDeg(), fitResult.LowerError(0) * TMath::RadToDeg(), mFitFunctor.laPreCorr[iDet] * TMath::RadToDeg(), vdFitResults[iDet], fitResult.LowerError(1), mFitFunctor.vdPreCorr[iDet], 100*fitResult.MinFcnValue());
170+
hVd->Fill(vdFitResults[iDet]);
171+
hLa->Fill(laFitResults[iDet]* TMath::RadToDeg());
131172
}
132173
printf("-------- Finished fits\n");
174+
175+
std::cout<<"number of chambers with enough entries: "<<hVd->GetEntries()<<std::endl;;
176+
std::cout<<"vdrift mean: "<<hVd->GetMean()<<" sigma: "<<hVd->GetStdDev()<<std::endl;
177+
std::cout<<"lorentz angle mean: "<<hLa->GetMean()<<" sigma: "<<hLa->GetStdDev()<<std::endl;
178+
179+
TCanvas* c = new TCanvas("c", "c", 1200, 700);
180+
c->Divide(2,1);
181+
c->cd(1);
182+
hVd->SetLineColor(kBlue); hVd->SetLineWidth(2);
183+
hVd->Draw();
184+
c->cd(2);
185+
hLa->SetLineColor(kBlue); hLa->SetLineWidth(2);
186+
hLa->Draw();
187+
c->SaveAs("VdExBValues.pdf");
188+
133189

134190
//----------------------------------------------------
135191
// Write

0 commit comments

Comments
 (0)