Project
Loading...
Searching...
No Matches
MatLayerCylSet.cxx
Go to the documentation of this file.
1// Copyright 2019-2020 CERN and copyright holders of ALICE O2.
2// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders.
3// All rights not expressly granted are reserved.
4//
5// This software is distributed under the terms of the GNU General Public
6// License v3 (GPL Version 3), copied verbatim in the file "COPYING".
7//
8// In applying this license CERN does not waive the privileges and immunities
9// granted to it by virtue of its status as an Intergovernmental Organization
10// or submit itself to any jurisdiction.
13
16
17#ifndef GPUCA_ALIGPUCODE // this part is unvisible on GPU version
18#include "GPUCommonLogger.h"
19#include <TFile.h>
21//#define _DBG_LOC_ // for local debugging only
22
23#endif // !GPUCA_ALIGPUCODE
24#undef NDEBUG
25using namespace o2::base;
26
28
29#ifndef GPUCA_ALIGPUCODE // this part is unvisible on GPU version
30
31//________________________________________________________________________________
32void MatLayerCylSet::addLayer(float rmin, float rmax, float zmax, float dz, float drphi)
33{
34 // add new layer checking for overlaps
36 assert(rmin < rmax && zmax > 0 && dz > 0 && drphi > 0);
38 int nlr = getNLayers();
39 if (!nlr) {
40 // book local storage
41 auto sz = sizeof(MatLayerCylSetLayout);
44 mFlatBufferSize = sz;
45 //--------------????
46 get()->mRMin = 1.e99;
47 get()->mRMax = 0.;
48 }
49
50 for (int il = 0; il < nlr; il++) {
51 const auto& lr = getLayer(il);
52 if (lr.getRMax() > rmin && rmax > lr.getRMin()) {
53 LOG(fatal) << "new layer overlaps with layer " << il;
54 }
55 }
56 auto* oldLayers = o2::gpu::resizeArray(get()->mLayers, nlr, nlr + 1);
57 // dynamyc buffers of old layers were used in new ones, detach them
58 for (int i = nlr; i--;) {
59 oldLayers[i].clearInternalBufferPtr();
60 }
61 delete[] oldLayers;
62 get()->mLayers[nlr].initSegmentation(rmin, rmax, zmax, dz, drphi);
63 get()->mNLayers++;
64 get()->mRMin = get()->mRMin > rmin ? rmin : get()->mRMin;
65 get()->mRMax = get()->mRMax < rmax ? rmax : get()->mRMax;
66 get()->mZMax = get()->mZMax < zmax ? zmax : get()->mZMax;
67 get()->mRMin2 = get()->mRMin * get()->mRMin;
68 get()->mRMax2 = get()->mRMax * get()->mRMax;
69}
70
71//________________________________________________________________________________
73{
76
77 int nlr = getNLayers();
78 if (!nlr) {
79 LOG(error) << "The LUT is not yet initialized";
80 return;
81 }
82 if (get()->mR2Intervals) {
83 LOG(error) << "The LUT is already populated";
84 return;
85 }
86 for (int i = 0; i < nlr; i++) {
87 printf("Populating with %d trials Lr %3d ", ntrPerCell, i);
88 get()->mLayers[i].print();
89 get()->mLayers[i].populateFromTGeo(ntrPerCell);
90 }
92}
93
94//________________________________________________________________________________
96{
97 // build layer search structures
99 int nlr = getNLayers();
100 int nR2Int = 2 * (nlr + 1);
101 o2::gpu::resizeArray(get()->mR2Intervals, 0, nR2Int);
102 o2::gpu::resizeArray(get()->mInterval2LrID, 0, nR2Int);
103 get()->mR2Intervals[0] = get()->mRMin2;
104 get()->mR2Intervals[1] = get()->mRMax2;
105 get()->mInterval2LrID[0] = 0;
106 auto& nRIntervals = get()->mNRIntervals;
107 nRIntervals = 1;
108
109 for (int i = 1; i < nlr; i++) {
110 const auto& lr = getLayer(i);
111 if (o2::math_utils::sqrt(lr.getRMin2()) > o2::math_utils::sqrt(get()->mR2Intervals[nRIntervals] + Ray::Tiny)) {
112 // register gap
113 get()->mInterval2LrID[nRIntervals] = -1;
114 get()->mR2Intervals[++nRIntervals] = lr.getRMin2();
115 }
116 get()->mInterval2LrID[nRIntervals] = i;
117 get()->mR2Intervals[++nRIntervals] = lr.getRMax2();
118 }
119 delete[] o2::gpu::resizeArray(get()->mInterval2LrID, nR2Int, nRIntervals); // rebook with precise size
120 delete[] o2::gpu::resizeArray(get()->mR2Intervals, nR2Int, ++nRIntervals); // rebook with precise size
121 //
122}
123
124//________________________________________________________________________________
125void MatLayerCylSet::dumpToTree(const std::string& outName) const
126{
128
129 o2::utils::TreeStreamRedirector dump(outName.data(), "recreate");
130 for (int i = 0; i < getNLayers(); i++) {
131 const auto& lr = getLayer(i);
132 float r = 0.5 * (lr.getRMin() + lr.getRMax());
133 // per cell dump
134 int nphib = lr.getNPhiBins();
135 for (int ip = 0; ip < nphib; ip++) {
136 float phi = 0.5 * (lr.getPhiBinMin(ip) + lr.getPhiBinMax(ip));
137 float sn, cs;
138 int ips = lr.phiBin2Slice(ip);
139 char merge = 0; // not mergeable
140 if (ip + 1 < nphib) {
141 int ips1 = lr.phiBin2Slice(ip + 1);
142 merge = ips == ips1 ? -1 : lr.canMergePhiSlices(ips, ips1); // -1 for already merged
143 } else {
144 merge = -2; // last one
145 }
146 o2::math_utils::sincos(phi, sn, cs);
147 float x = r * cs, y = r * sn;
148 for (int iz = 0; iz < lr.getNZBins(); iz++) {
149 float z = 0.5 * (lr.getZBinMin(iz) + lr.getZBinMax(iz));
150 auto cell = lr.getCellPhiBin(ip, iz);
151 dump << "cell"
152 << "ilr=" << i << "r=" << r << "phi=" << phi << "x=" << x << "y=" << y << "z=" << z << "ip=" << ip << "ips=" << ips << "iz=" << iz
153 << "mrgnxt=" << merge << "val=" << cell << "\n";
154 }
155 }
156 //
157 // statistics per layer
158 MatCell mean, rms;
159 lr.getMeanRMS(mean, rms);
160 dump << "lay"
161 << "ilr=" << i << "r=" << r << "mean=" << mean << "rms=" << rms << "\n";
162 }
163}
164
165//________________________________________________________________________________
166void MatLayerCylSet::writeToFile(const std::string& outFName)
167{
169
170 TFile outf(outFName.data(), "recreate");
171 if (outf.IsZombie()) {
172 return;
173 }
174 outf.WriteObjectAny(this, Class(), "ccdb_object");
175 outf.Close();
176}
177
179{
181 LOG(info) << "Layer voxel already initialized; Aborting";
182 return;
183 }
184 LOG(info) << "Initializing voxel layer lookup";
185 // do some check if voxels are dimensioned correctly
186 if (LayerRMax < get()->mRMax) {
187 LOG(fatal) << "Cannot initialized layer voxel lookup due to dimension problem (fix constants in MatLayerCylSet.h)";
188 }
189 for (int voxel = 0; voxel < NumVoxels; ++voxel) {
190 // check the 2 extremes of this voxel "covering"
191 const auto lowerR = voxel * VoxelRDelta;
192 const auto upperR = lowerR + VoxelRDelta;
193 const auto lowerSegment = searchSegment(lowerR * lowerR);
194 const auto upperSegment = searchSegment(upperR * upperR);
195 mLayerVoxelLU[2 * voxel] = lowerSegment;
196 mLayerVoxelLU[2 * voxel + 1] = upperSegment;
197 }
199}
200
201//________________________________________________________________________________
202MatLayerCylSet* MatLayerCylSet::loadFromFile(const std::string& inpFName)
203{
204 TFile inpf(inpFName.data());
205 if (inpf.IsZombie()) {
206 LOG(error) << "Failed to open input file " << inpFName;
207 return nullptr;
208 }
209 MatLayerCylSet* mb = reinterpret_cast<MatLayerCylSet*>(inpf.GetObjectChecked("ccdb_object", Class()));
210 if (!mb && !(mb = reinterpret_cast<MatLayerCylSet*>(inpf.GetObjectChecked("MatBud", Class())))) { // for old objects
211 LOG(error) << "Failed to load mat.LUT from " << inpFName;
212 return nullptr;
213 }
214 auto rptr = rectifyPtrFromFile(mb);
215 return rptr;
216}
217
218//________________________________________________________________________________
220{
221 // rectify object loaded from file
222 if (ptr && !ptr->get()) {
223 ptr->fixPointers();
224 }
225 ptr->initLayerVoxelLU();
226 return ptr;
227}
228
229//________________________________________________________________________________
231{
232 // merge similar (whose relative budget does not differ within maxRelDiff) phi slices
233 assert(mConstructionMask == InProgress);
234 for (int i = getNLayers(); i--;) {
235 get()->mLayers[i].optimizePhiSlices(maxRelDiff);
236 }
237 // flatten(); // RS: TODO
238}
239
240//________________________________________________________________________________
242{
244 if (!get()) {
245 printf("Not initialized yet\n");
246 return;
247 }
249 LOG(warning) << "Object is not yet flattened";
250 }
251 for (int i = 0; i < getNLayers(); i++) {
252 printf("#%3d | ", i);
254 }
255 printf("%.2f < R < %.2f %d layers with total size %.2f MB\n", getRMin(), getRMax(), getNLayers(),
256 float(getFlatBufferSize()) / 1024 / 1024);
257}
258
259//________________________________________________________________________________
260void MatLayerCylSet::scaleLayersByID(int lrFrom, int lrTo, float factor, bool _x2x0, bool _rho)
261{
262 lrFrom = std::max(0, std::min(lrFrom, get()->mNLayers - 1));
263 lrTo = std::max(0, std::min(lrTo, get()->mNLayers - 1));
264 int dir = lrFrom >= lrTo ? -1 : 1;
265 lrTo += dir;
266 for (int i = lrFrom; i != lrTo; i += dir) {
267 get()->mLayers[i].scale(factor, _x2x0, _rho);
268 }
269}
270
271//________________________________________________________________________________
272void MatLayerCylSet::scaleLayersByR(float rFrom, float rTo, float factor, bool _x2x0, bool _rho)
273{
274 if (rFrom > rTo) {
275 std::swap(rFrom, rTo);
276 }
277 Ray ray(std::max(getRMin(), rFrom), 0., 0., std::min(getRMax(), rTo), 0., 0.);
278 short lmin, lmax;
279 if (!getLayersRange(ray, lmin, lmax)) {
280 LOGP(warn, "No layers found for {} < r < {}", rFrom, rTo);
281 return;
282 }
283 scaleLayersByID(lmin, lmax, factor, _x2x0, _rho);
284}
285
286#endif
287
288#ifndef GPUCA_GPUCODE
289//________________________________________________________________________________
291{
292 std::size_t sz = alignSize(sizeof(MatLayerCylSetLayout), getBufferAlignmentBytes()); // hold data members
293
294 sz = alignSize(sz + get()->mNLayers * sizeof(MatLayerCyl), MatLayerCyl::getClassAlignmentBytes());
295 sz = alignSize(sz + (get()->mNRIntervals + 1) * sizeof(float), getBufferAlignmentBytes());
296 sz = alignSize(sz + get()->mNRIntervals * sizeof(int), getBufferAlignmentBytes());
297
298 for (int i = 0; i < getNLayers(); i++) {
300 }
301 return sz;
302}
303#endif // ! GPUCA_GPUCODE
304
305//_________________________________________________________________________________________________
306GPUd() MatBudget MatLayerCylSet::getMatBudget(float x0, float y0, float z0, float x1, float y1, float z1) const
307{
308 // get material budget traversed on the line between point0 and point1
309 MatBudget rval;
310 Ray ray(x0, y0, z0, x1, y1, z1);
311 short lmin, lmax; // get innermost and outermost relevant layer
312 if (ray.isTooShort() || !getLayersRange(ray, lmin, lmax)) {
313 rval.length = ray.getDist();
314 return rval;
315 }
316 short lrID = lmax;
317 while (lrID >= lmin) { // go from outside to inside
318 const auto& lr = getLayer(lrID);
319 int nphiSlices = lr.getNPhiSlices();
320 int nc = ray.crossLayer(lr); // determines how many crossings this ray has with this tubular layer
321 for (int ic = nc; ic--;) {
322 float cross1, cross2;
323 ray.getCrossParams(ic, cross1, cross2); // tmax,tmin of crossing the layer
324
325 auto phi0 = ray.getPhi(cross1), phi1 = ray.getPhi(cross2), dPhi = phi0 - phi1;
326 auto phiID = lr.getPhiSliceID(phi0), phiIDLast = lr.getPhiSliceID(phi1);
327 // account for eventual wrapping around 0
328 if (dPhi > 0.f) {
329 if (dPhi > o2::constants::math::PI) { // wraps around phi=0
330 phiIDLast += nphiSlices;
331 }
332 } else {
333 if (dPhi < -o2::constants::math::PI) { // wraps around phi=0
334 phiID += nphiSlices;
335 }
336 }
337
338 int stepPhiID = phiID > phiIDLast ? -1 : 1;
339 bool checkMorePhi = true;
340 auto tStartPhi = cross1, tEndPhi = 0.f;
341 do {
342 // get the path in the current phi slice
343 if (phiID == phiIDLast) {
344 tEndPhi = cross2;
345 checkMorePhi = false;
346 } else { // last phi slice still not reached
347 tEndPhi = ray.crossRadial(lr, (stepPhiID > 0 ? phiID + 1 : phiID) % nphiSlices);
348 if (tEndPhi == Ray::InvalidT) {
349 break; // ray parallel to radial line, abandon check for phi bin change
350 }
351 }
352 auto zID = lr.getZBinID(ray.getZ(tStartPhi));
353 auto zIDLast = lr.getZBinID(ray.getZ(tEndPhi));
354 // check if Zbins are crossed
355
356#ifdef _DBG_LOC_
357 printf("-- Zdiff (%3d : %3d) mode: t: %+e %+e\n", zID, zIDLast, tStartPhi, tEndPhi);
358#endif
359
360 if (zID != zIDLast) {
361 auto stepZID = zID < zIDLast ? 1 : -1;
362 bool checkMoreZ = true;
363 auto tStartZ = tStartPhi, tEndZ = 0.f;
364 do {
365 if (zID == zIDLast) {
366 tEndZ = tEndPhi;
367 checkMoreZ = false;
368 } else {
369 tEndZ = ray.crossZ(lr.getZBinMin(stepZID > 0 ? zID + 1 : zID));
370 if (tEndZ == Ray::InvalidT) { // track normal to Z axis, abandon Zbin change test
371 break;
372 }
373 }
374 // account materials of this step
375 float step = tEndZ > tStartZ ? tEndZ - tStartZ : tStartZ - tEndZ; // the real step is ray.getDist(tEnd-tStart), will rescale all later
376 const auto& cell = lr.getCell(phiID % nphiSlices, zID);
377 rval.meanRho += cell.meanRho * step;
378 rval.meanX2X0 += cell.meanX2X0 * step;
379 rval.length += step;
380
381#ifdef _DBG_LOC_
382 float pos0[3] = {ray.getPos(tStartZ, 0), ray.getPos(tStartZ, 1), ray.getPos(tStartZ, 2)};
383 float pos1[3] = {ray.getPos(tEndZ, 0), ray.getPos(tEndZ, 1), ray.getPos(tEndZ, 2)};
384 printf(
385 "Lr#%3d / cross#%d : account %f<t<%f at phiSlice %d | Zbin: %3d (%3d) |[%+e %+e +%e]:[%+e %+e %+e] "
386 "Step: %.3e StrpCor: %.3e\n",
387 lrID, ic, tEndZ, tStartZ, phiID % nphiSlices, zID, zIDLast,
388 pos0[0], pos0[1], pos0[2], pos1[0], pos1[1], pos1[2], step, ray.getDist(step));
389#endif
390
391 tStartZ = tEndZ;
392 zID += stepZID;
393 } while (checkMoreZ);
394 } else {
395 float step = tEndPhi > tStartPhi ? tEndPhi - tStartPhi : tStartPhi - tEndPhi; // the real step is |ray.getDist(tEnd-tStart)|, will rescale all later
396 const auto& cell = lr.getCell(phiID % nphiSlices, zID);
397 rval.meanRho += cell.meanRho * step;
398 rval.meanX2X0 += cell.meanX2X0 * step;
399 rval.length += step;
400
401#ifdef _DBG_LOC_
402 float pos0[3] = {ray.getPos(tStartPhi, 0), ray.getPos(tStartPhi, 1), ray.getPos(tStartPhi, 2)};
403 float pos1[3] = {ray.getPos(tEndPhi, 0), ray.getPos(tEndPhi, 1), ray.getPos(tEndPhi, 2)};
404 printf(
405 "Lr#%3d / cross#%d : account %f<t<%f at phiSlice %d | Zbin: %3d ----- |[%+e %+e +%e]:[%+e %+e %+e]"
406 "Step: %.3e StrpCor: %.3e\n",
407 lrID, ic, tEndPhi, tStartPhi, phiID % nphiSlices, zID,
408 pos0[0], pos0[1], pos0[2], pos1[0], pos1[1], pos1[2], step, ray.getDist(step));
409#endif
410 }
411 //
412 tStartPhi = tEndPhi;
413 phiID += stepPhiID;
414
415 } while (checkMorePhi);
416 }
417 lrID--;
418 } // loop over layers
419
420 if (rval.length != 0.f) {
421 rval.meanRho /= rval.length; // average
422 rval.meanX2X0 *= ray.getDist(); // normalize
423 }
424 rval.length = ray.getDist();
425
426#ifdef _DBG_LOC_
427 printf("<rho> = %e, x2X0 = %e | step = %e\n", rval.meanRho, rval.meanX2X0, rval.length);
428#endif
429 return rval;
430}
431
432//_________________________________________________________________________________________________
433GPUd() bool MatLayerCylSet::getLayersRange(const Ray& ray, short& lmin, short& lmax) const
434{
435 // get range of layers corresponding to rmin/rmax
436 //
437 lmin = lmax = -1;
438 float rmin2, rmax2;
439 ray.getMinMaxR2(rmin2, rmax2);
440
441 if (rmin2 >= getRMax2() || rmax2 <= getRMin2()) {
442 return false;
443 }
444 int lmxInt, lmnInt;
445 if (!mInitializedLayerVoxelLU) {
446 lmxInt = rmax2 < getRMax2() ? searchSegment(rmax2, 0) : get()->mNRIntervals - 2;
447 lmnInt = rmin2 >= getRMin2() ? searchSegment(rmin2, 0, lmxInt + 1) : 0;
448 } else {
449 lmxInt = rmax2 < getRMax2() ? searchLayerFast(rmax2, 0) : get()->mNRIntervals - 2;
450 lmnInt = rmin2 >= getRMin2() ? searchLayerFast(rmin2, 0, lmxInt + 1) : 0;
451 }
452
453 const auto* interval2LrID = get()->mInterval2LrID;
454 lmax = interval2LrID[lmxInt];
455 lmin = interval2LrID[lmnInt];
456 // make sure lmnInt and/or lmxInt are not in the gap
457 if (lmax < 0) {
458 lmax = interval2LrID[lmxInt - 1]; // rmax2 is in the gap, take highest layer below rmax2
459 }
460 if (lmin < 0) {
461 lmin = interval2LrID[lmnInt + 1]; // rmin2 is in the gap, take lowest layer above rmin2
462 }
463 return lmin <= lmax; // valid if both are not in the same gap
464}
465
466GPUd() int MatLayerCylSet::searchLayerFast(float r2, int low, int high) const
467{
468 // we can avoid the sqrt .. at the cost of more memory in the lookup
469 const auto index = 2 * int(o2::gpu::CAMath::Sqrt(r2) * InvVoxelRDelta);
470 const auto layersfirst = mLayerVoxelLU[index];
471 const auto layerslast = mLayerVoxelLU[index + 1];
472 if (layersfirst != layerslast) {
473 // this means the voxel is undecided and we revert to search
474 return searchSegment(r2, layersfirst, layerslast + 1);
475 }
476 return layersfirst;
477}
478
479GPUd() int MatLayerCylSet::searchSegment(float val, int low, int high) const
480{
482 if (low < 0) {
483 low = 0;
484 }
485 if (high < 0) {
486 high = get()->mNRIntervals;
487 }
488 int mid = (low + high) >> 1;
489 const auto* r2Intervals = get()->mR2Intervals;
490 while (mid != low) {
491 if (val < r2Intervals[mid]) {
492 high = mid;
493 } else {
494 low = mid;
495 }
496 mid = (low + high) >> 1;
497 }
498
499 return mid;
500}
501
502#ifndef GPUCA_ALIGPUCODE // this part is unvisible on GPU version
503
505{
506 // make object flat: move all content to single internally allocated buffer
507 assert(mConstructionMask == InProgress);
508
509 int sz = estimateFlatBufferSize();
510 // create new internal buffer with total size and copy data
513 mFlatBufferSize = sz;
514 int nLr = getNLayers();
515
516 auto offs = alignSize(sizeof(MatLayerCylSetLayout), getBufferAlignmentBytes()); // account for the alignment
517 // move array of layer pointers to the flat array
518 auto* oldLayers = o2::gpu::resizeArray(get()->mLayers, nLr, nLr, (MatLayerCyl*)(mFlatBufferPtr + offs));
519 // dynamyc buffers of old layers were used in new ones, detach them
520 for (int i = nLr; i--;) {
521 oldLayers[i].clearInternalBufferPtr();
522 }
523 delete[] oldLayers;
524 offs = alignSize(offs + nLr * sizeof(MatLayerCyl), MatLayerCyl::getClassAlignmentBytes()); // account for the alignment
525
526 // move array of R2 boundaries to the flat array
527 delete[] o2::gpu::resizeArray(get()->mR2Intervals, nLr + 1, nLr + 1, (float*)(mFlatBufferPtr + offs));
528 offs = alignSize(offs + (nLr + 1) * sizeof(float), getBufferAlignmentBytes()); // account for the alignment
529
530 // move array of R2 boundaries to the flat array
531 delete[] o2::gpu::resizeArray(get()->mInterval2LrID, nLr, nLr, (int*)(mFlatBufferPtr + offs));
532 offs = alignSize(offs + nLr * sizeof(int), getBufferAlignmentBytes()); // account for the alignment
533
534 for (int il = 0; il < nLr; il++) {
535 MatLayerCyl& lr = get()->mLayers[il];
536 lr.flatten(mFlatBufferPtr + offs);
537 offs = alignSize(offs + lr.getFlatBufferSize(), getBufferAlignmentBytes()); // account for the alignment
538 }
540}
541
542//______________________________________________
543void MatLayerCylSet::moveBufferTo(char* newFlatBufferPtr)
544{
546 flatObject::moveBufferTo(newFlatBufferPtr);
548}
549#endif // !GPUCA_ALIGPUCODE
550
551#ifndef GPUCA_GPUCODE
552//______________________________________________
553void MatLayerCylSet::setFutureBufferAddress(char* futureFlatBufferPtr)
554{
557 fixPointers(mFlatBufferPtr, futureFlatBufferPtr, false); // flag that futureFlatBufferPtr is not valid yet
558 flatObject::setFutureBufferAddress(futureFlatBufferPtr);
559}
560
561//______________________________________________
562void MatLayerCylSet::setActualBufferAddress(char* actualFlatBufferPtr)
563{
566 fixPointers(actualFlatBufferPtr);
567}
568//______________________________________________
569void MatLayerCylSet::cloneFromObject(const MatLayerCylSet& obj, char* newFlatBufferPtr)
570{
572 flatObject::cloneFromObject(obj, newFlatBufferPtr);
574}
575
576//______________________________________________
577void MatLayerCylSet::fixPointers(char* newBasePtr)
578{
579 // fix pointers on the internal structure of the flat buffer after retrieving it from the file
580 if (newBasePtr) {
581 mFlatBufferPtr = newBasePtr; // used to impose external pointer
582 } else {
583 mFlatBufferPtr = mFlatBufferContainer; // impose pointer after reading from file
584 }
585 auto offs = alignSize(sizeof(MatLayerCylSetLayout), getBufferAlignmentBytes()); // account for the alignment
586 char* newPtr = mFlatBufferPtr + offs; // correct pointer on MatLayerCyl*
587 char* oldPtr = reinterpret_cast<char*>(get()->mLayers); // old pointer read from the file
588 fixPointers(oldPtr, newPtr);
589}
590
591//______________________________________________
592void MatLayerCylSet::fixPointers(char* oldPtr, char* newPtr, bool newPtrValid)
593{
594 // fix pointers on the internal structure of the flat buffer after retrieving it from the file
595 auto* layPtr = get()->mLayers;
596 get()->mLayers = flatObject::relocatePointer(oldPtr, newPtr, get()->mLayers);
597 get()->mR2Intervals = flatObject::relocatePointer(oldPtr, newPtr, get()->mR2Intervals);
598 get()->mInterval2LrID = flatObject::relocatePointer(oldPtr, newPtr, get()->mInterval2LrID);
599 if (newPtrValid) {
600 layPtr = get()->mLayers;
601 }
602 for (int i = 0; i < getNLayers(); i++) {
603 layPtr[i].setFlatPointer(flatObject::relocatePointer(oldPtr, newPtr, layPtr[i].getFlatBufferPtr()));
604 layPtr[i].fixPointers(oldPtr, newPtr);
605 }
606}
607#endif // !GPUCA_GPUCODE
608
609#ifndef GPUCA_ALIGPUCODE // this part is unvisible on GPU version
610
611MatLayerCylSet* MatLayerCylSet::extractCopy(float rmin, float rmax, float tolerance) const
612{
613 Ray ray(std::max(getRMin(), rmin), 0., 0., std::min(getRMax(), rmax), 0., 0.);
614 short lmin, lmax;
615 if (!getLayersRange(ray, lmin, lmax)) {
616 LOGP(warn, "No layers found for {} < r < {}", rmin, rmax);
617 return nullptr;
618 }
619 LOGP(info, "Will extract layers {}:{} (out of {} layers) for {} < r < {}", lmin, lmax, getNLayers(), rmin, rmax);
620 MatLayerCylSet* copy = new MatLayerCylSet();
621 int lrCount = 0;
622 for (int il = lmin; il <= lmax; il++) {
623 const auto& lr = getLayer(il);
624 float drphi = lr.getDPhi() * (lr.getRMin() + lr.getRMax()) / 2. * 0.999;
625 copy->addLayer(lr.getRMin(), lr.getRMax(), lr.getZMax(), lr.getDZ(), drphi);
626 auto& lrNew = copy->getLayer(lrCount);
627 for (int iz = 0; iz < lrNew.getNZBins(); iz++) {
628 for (int ip = 0; ip < lrNew.getNPhiBins(); ip++) {
629 lrNew.getCellPhiBin(ip, iz).set(lr.getCellPhiBin(ip, iz));
630 }
631 }
632 lrCount++;
633 }
634
635 copy->finalizeStructures();
636 copy->optimizePhiSlices(tolerance);
637 copy->flatten();
638 return copy;
639}
640
641#endif
int32_t i
#define GPUd()
Declarations for the wrapper for the set of cylindrical material layers.
useful math constants
TBranch * ptr
void merge(Options const &options)
void setActualBufferAddress(char *actualFlatBufferPtr)
static constexpr size_t getBufferAlignmentBytes()
Gives minimal alignment in bytes required for the flat buffer.
void addLayer(float rmin, float rmax, float zmax, float dz, float drphi)
void optimizePhiSlices(float maxRelDiff=0.05)
void cloneFromObject(const MatLayerCylSet &obj, char *newFlatBufferPtr)
MatLayerCyl & getLayer(int i)
static MatLayerCylSet * loadFromFile(const std::string &inpFName="matbud.root")
static constexpr float VoxelRDelta
uint16_t mLayerVoxelLU[2 *NumVoxels]
static constexpr int NumVoxels
void print(bool data=false) const
void populateFromTGeo(int ntrPerCel=10)
bool mInitializedLayerVoxelLU
helper structure to lookup a layer based on known radius (static dimension for easy copy to GPU)
void scaleLayersByR(float rFrom, float rTo, float factor, bool _x2x0=true, bool _rho=true)
void moveBufferTo(char *newFlatBufferPtr)
MatLayerCylSet * extractCopy(float rmin, float rmax, float tol=1e-3) const
std::size_t estimateFlatBufferSize() const
GPUCA_ALIGPUCODE.
void dumpToTree(const std::string &outName="matbudTree.root") const
static MatLayerCylSet * rectifyPtrFromFile(MatLayerCylSet *ptr)
void setFutureBufferAddress(char *futureFlatBufferPtr)
static constexpr float LayerRMax
void scaleLayersByID(int lrFrom, int lrTo, float factor, bool _x2x0=true, bool _rho=true)
void fixPointers(char *newPtr=nullptr)
void writeToFile(const std::string &outFName="matbud.root")
static constexpr size_t getClassAlignmentBytes()
Gives minimal alignment in bytes required for the class object.
void print(bool data=false) const
void flatten(char *newPtr)
MatCell & getCellPhiBin(int iphi, int iz)
Definition MatLayerCyl.h:96
static constexpr float InvalidT
Definition Ray.h:46
static constexpr float Tiny
Definition Ray.h:47
GPUCA_GPUCODE.
Definition FlatObject.h:176
void setFutureBufferAddress(char *futureFlatBufferPtr)
Definition FlatObject.h:557
uint32_t mConstructionMask
mask for constructed object members, first two bytes are used by this class
Definition FlatObject.h:323
int32_t mFlatBufferSize
size of the flat buffer
Definition FlatObject.h:322
char * mFlatBufferContainer
Definition FlatObject.h:324
static size_t alignSize(size_t sizeBytes, size_t alignmentBytes)
_______________ Generic utilities _______________________________________________
Definition FlatObject.h:277
static T * relocatePointer(const char *oldBase, char *newBase, const T *ptr)
Relocates a pointer inside a buffer to the new buffer address.
Definition FlatObject.h:285
void moveBufferTo(char *newBufferPtr)
Definition FlatObject.h:396
void cloneFromObject(const FlatObject &obj, char *newFlatBufferPtr)
Definition FlatObject.h:373
size_t getFlatBufferSize() const
Gives size of the flat buffer.
Definition FlatObject.h:256
@ InProgress
construction started: temporary memory is reserved
Definition FlatObject.h:319
@ Constructed
the object is constructed, temporary memory is released
Definition FlatObject.h:318
const char * getFlatBufferPtr() const
Gives pointer to the flat buffer.
Definition FlatObject.h:259
void dump(const std::string what, DPMAP m, int verbose)
Definition dcs-ccdb.cxx:79
GLint GLenum GLint x
Definition glcorearb.h:403
GLuint GLfloat GLfloat GLfloat GLfloat y1
Definition glcorearb.h:5034
GLuint GLfloat GLfloat GLfloat x1
Definition glcorearb.h:5034
GLuint index
Definition glcorearb.h:781
GLboolean * data
Definition glcorearb.h:298
GLuint GLfloat x0
Definition glcorearb.h:5034
GLuint GLfloat * val
Definition glcorearb.h:1582
GLboolean r
Definition glcorearb.h:1233
GLuint GLfloat GLfloat y0
Definition glcorearb.h:5034
GLdouble GLdouble GLdouble z
Definition glcorearb.h:843
constexpr float PI
T * resizeArray(T *&ptr, int32_t oldSize, int32_t newSize, T *newPtr=nullptr)
Definition FlatObject.h:121
auto get(const std::byte *buffer, size_t=0)
Definition DataHeader.h:454
float length
length in material
Definition MatCell.h:55
float meanRho
mean density, g/cm^3
Definition MatCell.h:30
float meanX2X0
fraction of radiaton lenght
Definition MatCell.h:31
LOG(info)<< "Compressed in "<< sw.CpuTime()<< " s"