Skip to content

Commit 869a001

Browse files
matthias-kleinerdavidrohr
matthias-kleiner
authored and
davidrohr
committed
Space-charge: Add conversion of IDCs to space-charge density
other changes: - Add function to set the space-charge density from two histograms (A and C-side) - Add function to convert IDCs to the space-charge density - Add function to add charge from other object (e.g. Adding space-charge density from primary ionization to space-charge density from IDCs) - Speed up rebinning by directly filling the space-charge density createTPCSpaceChargeCorrection - Fix initialization of fastTransform - Add protection when dumping to TTree IDCs - fix sign for drawing local coordinates
1 parent d0558ad commit 869a001

File tree

6 files changed

+218
-53
lines changed

6 files changed

+218
-53
lines changed

‎Detectors/TPC/calibration/src/IDCDrawHelper.cxx‎

Lines changed: 6 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -48,12 +48,13 @@ void o2::tpc::IDCDrawHelper::drawSector(const IDCDraw& idc, const unsigned int s
4848
can->SetTopMargin(0.04f);
4949
poly->Draw("colz");
5050

51+
const float sign = (Sector(sector).side() == Side::A) ? 1 : -1;
5152
for (unsigned int region = startRegion; region < endRegion; ++region) {
5253
for (unsigned int irow = 0; irow < Mapper::ROWSPERREGION[region]; ++irow) {
5354
for (unsigned int ipad = 0; ipad < Mapper::PADSPERROW[region][irow]; ++ipad) {
5455
const auto padNum = Mapper::getGlobalPadNumber(irow, ipad, region);
55-
const auto coordinate = coords[padNum];
56-
const float yPos = -static_cast<float>(coordinate.yVals[0] + coordinate.yVals[2]) / 2; // local coordinate system is mirrored
56+
const auto coordinate = coords[padNum]; // start from y=-13 to +13
57+
const float yPos = sign * static_cast<float>(coordinate.yVals[0] + coordinate.yVals[2]) / 2; // local coordinate system is mirrored for C side
5758
const float xPos = static_cast<float>(coordinate.xVals[0] + coordinate.xVals[2]) / 2;
5859
poly->Fill(xPos, yPos, idc.getIDC(sector, region, irow, ipad));
5960
}
@@ -120,11 +121,11 @@ TH2Poly* o2::tpc::IDCDrawHelper::drawSide(const IDCDraw& idc, const o2::tpc::Sid
120121
for (unsigned int ipad = 0; ipad < Mapper::PADSPERROW[region][irow]; ++ipad) {
121122
const auto padNum = Mapper::getGlobalPadNumber(irow, ipad, region);
122123
const float angDeg = 10.f + sector * 20;
123-
auto coordinate = coords[padNum];
124-
coordinate.rotate(angDeg);
124+
auto coordinate = coords[padNum]; // start from y=-13 to +13
125+
coordinate.rotate(angDeg); // start from y=1.2 to +28..
125126
const float yPos = static_cast<float>(coordinate.yVals[0] + coordinate.yVals[1] + coordinate.yVals[2] + coordinate.yVals[3]) / 4;
126127
const float xPos = static_cast<float>(coordinate.xVals[0] + coordinate.xVals[1] + coordinate.xVals[2] + coordinate.xVals[3]) / 4;
127-
const auto padTmp = getPad(ipad, region, irow, side);
128+
const auto padTmp = getPad(ipad, region, irow, side); // IDCs are in pad coordinates. pad0 A side: y=1.2, pad0 C side: y=28.1
128129
poly->Fill(xPos, yPos, idc.getIDC(sector, region, irow, padTmp));
129130
}
130131
}

‎Detectors/TPC/reconstruction/macro/createTPCSpaceChargeCorrection.C‎

Lines changed: 22 additions & 10 deletions
Original file line numberDiff line numberDiff line change
@@ -122,7 +122,8 @@ void createTPCSpaceChargeCorrection(
122122
spaceCharge->setGlobalCorrectionsFromFile(scFile, Side::C);
123123
TPCFastSpaceChargeCorrectionHelper::instance()->setGlobalSpaceChargeCorrection(getGlobalSpaceChargeCorrection);
124124

125-
std::unique_ptr<TPCFastTransform> fastTransform(TPCFastTransformHelperO2::instance()->create(0));
125+
std::unique_ptr<TPCFastSpaceChargeCorrection> spCorrection = TPCFastSpaceChargeCorrectionHelper::instance()->create();
126+
std::unique_ptr<TPCFastTransform> fastTransform(TPCFastTransformHelperO2::instance()->create(0, *spCorrection));
126127

127128
fastTransform->writeToFile(outputFileName, "ccdb_object");
128129

@@ -182,7 +183,8 @@ void createTPCSpaceChargeCorrectionLinearCombination(
182183
}
183184

184185
TPCFastSpaceChargeCorrectionHelper::instance()->setGlobalSpaceChargeCorrection(getGlobalSpaceChargeCorrectionLinearCombination);
185-
std::unique_ptr<TPCFastTransform> fastTransform(TPCFastTransformHelperO2::instance()->create(0));
186+
std::unique_ptr<TPCFastSpaceChargeCorrection> spCorrection = TPCFastSpaceChargeCorrectionHelper::instance()->create();
187+
std::unique_ptr<TPCFastTransform> fastTransform(TPCFastTransformHelperO2::instance()->create(0, *spCorrection));
186188
fastTransform->writeToFile(outputFileName, "ccdb_object");
187189

188190
if (debug > 0) {
@@ -230,7 +232,9 @@ void createTPCSpaceChargeCorrectionAnalytical(
230232

231233
TPCFastSpaceChargeCorrectionHelper::instance()->setLocalSpaceChargeCorrection(getLocalSpaceChargeCorrection);
232234

233-
std::unique_ptr<TPCFastTransform> fastTransform(TPCFastTransformHelperO2::instance()->create(0));
235+
std::unique_ptr<TPCFastSpaceChargeCorrection> spCorrection = TPCFastSpaceChargeCorrectionHelper::instance()->create();
236+
std::unique_ptr<TPCFastTransform> fastTransform(TPCFastTransformHelperO2::instance()->create(0, *spCorrection));
237+
234238
fastTransform->writeToFile(outputFileName, "ccdb_object");
235239

236240
if (debug > 0) {
@@ -425,7 +429,9 @@ void debugInterpolation(utils::TreeStreamRedirector& pcstream,
425429
// the original correction
426430
double gdC[3] = {0, 0, 0};
427431
Side side = slice < geo.getNumberOfSlicesA() ? Side::A : Side::C;
428-
spaceCharge->getCorrections(gx, gy, gz, side, gdC[0], gdC[1], gdC[2]);
432+
if (spaceCharge) {
433+
spaceCharge->getCorrections(gx, gy, gz, side, gdC[0], gdC[1], gdC[2]);
434+
}
429435
float ldxC, ldyC, ldzC;
430436
geo.convGlobalToLocal(slice, gdC[0], gdC[1], gdC[2], ldxC, ldyC, ldzC);
431437

@@ -437,18 +443,24 @@ void debugInterpolation(utils::TreeStreamRedirector& pcstream,
437443

438444
float pointCyl[3] = {gz, r, phi};
439445
double efield[3] = {0.0, 0.0, 0.0};
440-
double charge = spaceCharge->getDensityCyl(pointCyl[0], pointCyl[1], pointCyl[2], side);
441-
double potential = spaceCharge->getPotentialCyl(pointCyl[0], pointCyl[1], pointCyl[2], side);
442-
spaceCharge->getElectricFieldsCyl(pointCyl[0], pointCyl[1], pointCyl[2], side, efield[0], efield[1], efield[2]);
443-
spaceCharge->getLocalDistortionsCyl(pointCyl[0], pointCyl[1], pointCyl[2], side, ldD[0], ldD[1], ldD[2]);
444-
spaceCharge->getDistortions(gx, gy, gz, side, gdD[0], gdD[1], gdD[2]);
446+
double charge = 0;
447+
double potential = 0;
448+
if (spaceCharge) {
449+
charge = spaceCharge->getDensityCyl(pointCyl[0], pointCyl[1], pointCyl[2], side);
450+
potential = spaceCharge->getPotentialCyl(pointCyl[0], pointCyl[1], pointCyl[2], side);
451+
spaceCharge->getElectricFieldsCyl(pointCyl[0], pointCyl[1], pointCyl[2], side, efield[0], efield[1], efield[2]);
452+
spaceCharge->getLocalDistortionsCyl(pointCyl[0], pointCyl[1], pointCyl[2], side, ldD[0], ldD[1], ldD[2]);
453+
spaceCharge->getDistortions(gx, gy, gz, side, gdD[0], gdD[1], gdD[2]);
454+
}
445455

446456
double gD[3] = {gx + gdD[0], gy + gdD[1], gz + gdD[2]};
447457
float rD = std::sqrt(gD[0] * gD[0] + gD[1] * gD[1]);
448458

449459
// correction for the distorted point
450460
double gdDC[3] = {0, 0, 0};
451-
spaceCharge->getCorrections(gD[0], gD[1], gD[2], side, gdDC[0], gdDC[1], gdDC[2]);
461+
if (spaceCharge) {
462+
spaceCharge->getCorrections(gD[0], gD[1], gD[2], side, gdDC[0], gdDC[1], gdDC[2]);
463+
}
452464

453465
// exB
454466
float ldxC_ExB = 0, ldyC_ExB = 0, ldzC_ExB = 0;

‎Detectors/TPC/spacecharge/include/TPCSpaceCharge/DataContainer3D.h‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -185,6 +185,7 @@ struct DataContainer3D {
185185

186186
/// operator overload
187187
DataContainer3D<DataT>& operator*=(const DataT value);
188+
DataContainer3D<DataT>& operator+=(const DataContainer3D<DataT>& other);
188189

189190
private:
190191
unsigned short mZVertices{}; ///< number of z vertices

‎Detectors/TPC/spacecharge/include/TPCSpaceCharge/SpaceCharge.h‎

Lines changed: 40 additions & 8 deletions
Original file line numberDiff line numberDiff line change
@@ -118,10 +118,35 @@ class SpaceCharge
118118
/// \param hisSCDensity3D histogram for the space charge density
119119
void fillChargeDensityFromHisto(const TH3& hisSCDensity3D);
120120

121+
/// step 0: set the space-charge density from two TH3 histograms containing the space-charge density for A and C side seperately
122+
/// \param hisSCDensity3D_A histogram for the space charge density for A-side
123+
/// \param hisSCDensity3D_C histogram for the space charge density for C-side
124+
void fillChargeDensityFromHisto(const TH3& hisSCDensity3D_A, const TH3& hisSCDensity3D_C);
125+
126+
/// step 0: set the space-charge density from two TH3 histograms containing the space charge density for A and C side separately which are stored in a ROOT file
127+
/// \param file path to root file containing the space-charge density
128+
/// \param nameA name of the space-charge density histogram for the A-side
129+
/// \param nameC name of the space-charge density histogram for the C-side
130+
void fillChargeDensityFromHisto(const char* file, const char* nameA, const char* nameC);
131+
121132
/// step 0: set the space charge density from std::vector<CalDet> containing the space charge density. Each entry in the object corresponds to one z slice
122-
/// \param calSCDensity3D histogram for the space charge density
133+
/// \param calSCDensity3D pad-by-pad CalDet for the space charge density
123134
void fillChargeDensityFromCalDet(const std::vector<CalDet<float>>& calSCDensity3D);
124135

136+
/// Convert the IDCs to the number of ions for the ion backflow (primary ionization is not considered)
137+
/// \param idcZero map containing the IDCs values which will be converted to the space-charge density
138+
/// \param mapIBF map contains the pad-by-pad IBF in %
139+
/// \param ionDriftTimeMS ion drift time in ms
140+
/// \param normToPadArea normalize IDCs to pad area (should always be true as the normalization is performed in IDCFactorization::calcIDCZero
141+
static void convertIDCsToCharge(std::vector<CalDet<float>>& idcZero, const CalDet<float>& mapIBF, const float ionDriftTimeMS = 200, const bool normToPadArea = true);
142+
143+
/// Convert the IDCs to the number of ions for the ion backflow (primary ionization is not considered)
144+
/// \param idcZero map containing the IDCs which will be converted to the space-charge density
145+
/// \param mapIBF map contains the pad-by-pad IBF in %
146+
/// \param ionDriftTimeMS ion drift time in ms
147+
/// \param normToPadArea normalize IDCs to pad area (should always be true as the normalization is performed in IDCFactorization::calcIDCZero
148+
void fillChargeFromIDCs(std::vector<CalDet<float>>& idcZero, const CalDet<float>& mapIBF, const float ionDriftTimeMS = 200, const bool normToPadArea = true);
149+
125150
/// step 0: set the charge (number of ions) from std::vector<CalDet> containing the charge. Each entry in the object corresponds to one z slice.
126151
/// Normalization to the space charge is also done automatically
127152
/// \param calCharge3D histogram for the charge
@@ -252,6 +277,10 @@ class SpaceCharge
252277
/// \param side side for which the space-charge density will be scaled
253278
void scaleChargeDensity(const DataT scalingFactor, const Side side) { mDensity[side] *= scalingFactor; }
254279

280+
/// add space charge density from other object (this.mDensity = this.mDensity + other.mDensity)
281+
/// \param otherSC other space-charge object, which charge will be added to current object
282+
void addChargeDensity(const SpaceCharge<DataT>& otherSC);
283+
255284
/// step 3: calculate the local distortions and corrections with an electric field
256285
/// \param type calculate local corrections or local distortions: type = o2::tpc::SpaceCharge<>::Type::Distortions or o2::tpc::SpaceCharge<>::Type::Corrections
257286
/// \param formulaStruct struct containing a method to evaluate the electric field Er, Ez, Ephi (analytical formula or by TriCubic interpolator)
@@ -854,6 +883,16 @@ class SpaceCharge
854883
/// \param nRPoints number of vertices of the output in r
855884
/// \param nPhiPoints number of vertices of the output in phi
856885
/// \param randomize randomize points
886+
///
887+
/// Adding other TTree as friend:
888+
/// TFile f1("file_1.root");
889+
/// TTree* tree1 = (TTree*)f1.Get("tree");
890+
/// tree1->BuildIndex("globalIndex")
891+
///
892+
/// TFile f2("file_2.root");
893+
/// TTree* tree2 = (TTree*)f2.Get("tree");
894+
/// tree2->BuildIndex("globalIndex");
895+
/// tree1->AddFriend(tree2, "t2");
857896
void dumpToTree(const char* outFileName, const Side side, const int nZPoints = 50, const int nRPoints = 150, const int nPhiPoints = 180, const bool randomize = false) const;
858897

859898
/// dump to tree evaluated on the pads for given sector
@@ -1001,13 +1040,6 @@ class SpaceCharge
10011040
/// set phi coordinate between min phi max phi
10021041
DataT regulatePhi(const DataT posPhi, const Side side) const { return mGrid3D[side].clampToGridCircular(posPhi, 2); }
10031042

1004-
/// rebin the input space charge density histogram to desired binning
1005-
/// \param hOrig original histogram
1006-
/// \param nBinsZNew number of z bins of rebinned histogram
1007-
/// \param nBinsRNew number of r bins of rebinned histogram
1008-
/// \param nBinsPhiNew number of phi bins of rebinned histogram
1009-
static TH3DataT rebinDensityHisto(const TH3& hOrig, const unsigned short nBinsZNew, const unsigned short nBinsRNew, const unsigned short nBinsPhiNew);
1010-
10111043
/// function to calculate the drift paths of the electron whose starting position is delivered. Electric fields must be set!
10121044
/// \param elePos global position of the start position of the electron
10131045
/// \param nSamplingPoints number of output points of the electron drift path

‎Detectors/TPC/spacecharge/src/DataContainer3D.cxx‎

Lines changed: 7 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -233,6 +233,13 @@ DataContainer3D<DataT>& DataContainer3D<DataT>::operator*=(const DataT value)
233233
return *this;
234234
}
235235

236+
template <typename DataT>
237+
DataContainer3D<DataT>& DataContainer3D<DataT>::operator+=(const DataContainer3D<DataT>& other)
238+
{
239+
std::transform(mData.begin(), mData.end(), other.mData.begin(), mData.begin(), std::plus<>());
240+
return *this;
241+
}
242+
236243
template <typename DataT>
237244
size_t DataContainer3D<DataT>::getIndexZ(size_t index, const int nz, const int nr, const int nphi)
238245
{

0 commit comments

Comments
 (0)