Logo ROOT  

Reference Guide
 
Loading...
Searching...
No Matches
TPainter3dAlgorithms.cxx
Go to the documentation of this file.
1// @(#)root/histpainter:$Id$
2// Author: Rene Brun, Evgueni Tcherniaev, Olivier Couet 12/12/94
3
4/*************************************************************************
5 * Copyright (C) 1995-2000, Rene Brun and Fons Rademakers. *
6 * All rights reserved. *
7 * *
8 * For the licensing terms see $ROOTSYS/LICENSE. *
9 * For the list of contributors see $ROOTSYS/README/CREDITS. *
10 *************************************************************************/
11
12/*! \class TPainter3dAlgorithms
13 \ingroup Histpainter
14 \brief The Legos and Surfaces painter class.
15
163D graphics representations package.
17
18This package was originally written by Evgueni Tcherniaev from IHEP/Protvino.
19
20The original Fortran implementation was adapted to HIGZ/PAW by Olivier Couet
21and Evgueni Tcherniaev.
22
23This class is a subset of the original system. It has been converted to a C++
24class by Rene Brun.
25*/
26
27#include <cstdlib>
28
29#include "TROOT.h"
30#include "TPainter3dAlgorithms.h"
31#include "TVirtualPad.h"
32#include "THistPainter.h"
33#include "TH1.h"
34#include "TF3.h"
35#include "TView.h"
36#include "Hoption.h"
37#include "Hparam.h"
38#include "TMath.h"
39#include "TStyle.h"
40#include "THLimitsFinder.h"
41#include "TColor.h"
42
43const Double_t kRad = TMath::ATan(1)*Double_t(4)/Double_t(180);
44const Double_t kFdel = 0.;
45const Double_t kDel = 0.0001;
46const Int_t kNiso = 4;
47const Int_t kNmaxp = kNiso*13;
48const Int_t kNmaxt = kNiso*12;
49const Int_t kLmax = 12;
50const Int_t kF3FillColor1 = 201;
51const Int_t kF3FillColor2 = 202;
52const Int_t kF3LineColor = 203;
53
54extern TH1 *gCurrentHist; //these 3 globals should be replaced by class members
55extern Hoption_t Hoption;
56extern Hparam_t Hparam;
57
58////////////////////////////////////////////////////////////////////////////////
59/// Lego default constructor
60
61TPainter3dAlgorithms::TPainter3dAlgorithms(): TAttLine(1,1,1), TAttFill(1,0)
62{
63 Int_t i;
64 fNaphi = 0;
65 fIfrast = 0;
66 fMesh = 1;
67 fColorTop = 1;
68 fColorBottom = 1;
69 fEdgeIdx = -1;
70 fNlevel = 0;
71 fSystem = kCARTESIAN;
72 fDrawFace = nullptr;
73 fLegoFunction = nullptr;
74 fSurfaceFunction = nullptr;
75
76 TList *stack = nullptr;
77 if (gCurrentHist) stack = gCurrentHist->GetPainter()->GetStack();
78 fNStack = stack ? stack->GetSize() : 0;
79 fColorMain.resize(fNStack+1);
80 fColorDark.resize(fNStack+1);
81 fEdgeColor.resize(fNStack+1);
82 fEdgeStyle.resize(fNStack+1);
83 fEdgeWidth.resize(fNStack+1);
84
85 for (i=0;i<fNStack;i++) { fColorMain[i] = 1; fColorDark[i] = 1; fEdgeColor[i] = 1; fEdgeStyle[i] = 1; fEdgeWidth[i] = 1; }
86 for (i=0;i<3;i++) { fRmin[i] = 0; fRmax[i] = 1; }
87 for (i=0;i<4;i++) { fYls[i] = 0; }
88
89 for (i=0;i<30;i++) { fJmask[i] = 0; }
90 for (i=0;i<200;i++) { fLevelLine[i] = 0; }
91 for (i=0;i<465;i++) { fMask[i] = 0; }
92 for (i=0;i<258;i++) { fColorLevel[i] = 0; }
93 for (i=0;i<1200;i++) { fPlines[i] = 0.; }
94 for (i=0;i<200;i++) { fT[i] = 0.; }
95 for (i=0;i<2*NumOfSlices;i++) { fU[i] = 0.; fD[i] = 0.; }
96 for (i=0;i<12;i++) { fVls[i] = 0.; }
97 for (i=0;i<257;i++) { fFunLevel[i] = 0.; }
98 for (i=0;i<8;i++) { fF8[i] = 0.; }
99
100 fLoff = 0;
101 fNT = 0;
102 fNcolor = 0;
103 fNlines = 0;
104 fNqs = 0;
105 fNxrast = 0;
106 fNyrast = 0;
107 fIc1 = 0;
108 fIc2 = 0;
109 fIc3 = 0;
110 fQA = 0.;
111 fQD = 0.;
112 fQS = 0.;
113 fX0 = 0.;
114 fYdl = 0.;
115 fXrast = 0.;
116 fYrast = 0.;
117 fFmin = 0.;
118 fFmax = 0.;
119 fDXrast = 0.;
120 fDYrast = 0.;
121 fDX = 0.;
122}
123
124////////////////////////////////////////////////////////////////////////////////
125/// Normal default constructor
126///
127/// rmin[3], rmax[3] are the limits of the lego object depending on
128/// the selected coordinate system
129
130TPainter3dAlgorithms::TPainter3dAlgorithms(Double_t *rmin, Double_t *rmax, Int_t system)
131 : TAttLine(1,1,1), TAttFill(1,0)
132{
133 Int_t i;
134 Double_t psi;
135
136 fNaphi = 0;
137 fIfrast = 0;
138 fMesh = 1;
139 fColorTop = 1;
140 fColorBottom = 1;
141 fEdgeIdx = -1;
142 fNlevel = 0;
143 fSystem = system;
144 if (system == kCARTESIAN || system == kPOLAR) psi = 0;
145 else psi = 90;
146 fDrawFace = nullptr;
147 fLegoFunction = nullptr;
148 fSurfaceFunction = nullptr;
149
150 TList *stack = gCurrentHist->GetPainter()->GetStack();
151 fNStack = stack ? stack->GetSize() : 0;
152
153 fColorMain.resize(fNStack+1);
154 fColorDark.resize(fNStack+1);
155 fEdgeColor.resize(fNStack+1);
156 fEdgeStyle.resize(fNStack+1);
157 fEdgeWidth.resize(fNStack+1);
158
159 for (i=0;i<fNStack;i++) { fColorMain[i] = 1; fColorDark[i] = 1; fEdgeColor[i] = 1; fEdgeStyle[i] = 1; fEdgeWidth[i] = 1; }
160 for (i=0;i<3;i++) { fRmin[i] = rmin[i]; fRmax[i] = rmax[i]; }
161 for (i=0;i<4;i++) { fYls[i] = 0; }
162
163 for (i=0;i<30;i++) { fJmask[i] = 0; }
164 for (i=0;i<200;i++) { fLevelLine[i] = 0; }
165 for (i=0;i<465;i++) { fMask[i] = 0; }
166 for (i=0;i<258;i++) { fColorLevel[i] = 0; }
167 for (i=0;i<1200;i++) { fPlines[i] = 0.; }
168 for (i=0;i<200;i++) { fT[i] = 0.; }
169 for (i=0;i<2*NumOfSlices;i++) { fU[i] = 0.; fD[i] = 0.; }
170 for (i=0;i<12;i++) { fVls[i] = 0.; }
171 for (i=0;i<257;i++) { fFunLevel[i] = 0.; }
172 for (i=0;i<8;i++) { fF8[i] = 0.; }
173
174 fLoff = 0;
175 fNT = 0;
176 fNcolor = 0;
177 fNlines = 0;
178 fNqs = 0;
179 fNxrast = 0;
180 fNyrast = 0;
181 fIc1 = 0;
182 fIc2 = 0;
183 fIc3 = 0;
184 fQA = 0.;
185 fQD = 0.;
186 fQS = 0.;
187 fX0 = 0.;
188 fYdl = 0.;
189 fXrast = 0.;
190 fYrast = 0.;
191 fFmin = 0.;
192 fFmax = 0.;
193 fDXrast = 0.;
194 fDYrast = 0.;
195 fDX = 0.;
196
197 TView *view = gPad ? gPad->GetView() : nullptr;
198 if (!view)
199 view = TView::CreateView(fSystem, rmin, rmax);
200 if (view) {
201 view->SetView(gPad->GetPhi(), gPad->GetTheta(), psi, i);
202 view->SetRange(rmin,rmax);
203 }
204}
205
206////////////////////////////////////////////////////////////////////////////////
207/// destructor
208
212
213////////////////////////////////////////////////////////////////////////////////
214/// Draw back surfaces of surrounding box
215///
216/// \param[in] ang angle between X and Y axis
217
218void TPainter3dAlgorithms::BackBox(Double_t ang)
219{
220 static Int_t iface1[4] = { 1, 4, 8, 5 };
221 static Int_t iface2[4] = { 4, 3, 7, 8 };
222
223 TView *view = gPad ? gPad->GetView() : nullptr;
224 if (!view) {
225 Error("BackBox", "no TView in current pad");
226 return;
227 }
228
229 // Get corners of surrounding box
230 Double_t r[3*8], av[3*8];
231 Int_t ix1, ix2, iy1, iy2, iz1, iz2;
232 Double_t cosa = TMath::Cos(kRad*ang);
233 Double_t sina = TMath::Sin(kRad*ang);
234 view->AxisVertex(ang, av, ix1, ix2, iy1, iy2, iz1, iz2);
235 for (Int_t i = 0; i < 8; ++i) {
236 r[i*3 + 0] = av[i*3 + 0] + av[i*3 + 1]*cosa;
237 r[i*3 + 1] = av[i*3 + 1]*sina;
238 r[i*3 + 2] = av[i*3 + 2];
239 }
240
241 // Draw back faces
242 Int_t icodes[3] = { 0, 0, 0 };
243 Double_t tt[4];
244 tt[0] = r[(iface1[0]-1)*3 + 2];
245 tt[1] = r[(iface1[1]-1)*3 + 2];
246 tt[2] = r[(iface1[2]-1)*3 + 2];
247 tt[3] = r[(iface1[3]-1)*3 + 2];
248 (this->*fDrawFace)(icodes, r, 4, iface1, tt);
249 tt[0] = r[(iface2[0]-1)*3 + 2];
250 tt[1] = r[(iface2[1]-1)*3 + 2];
251 tt[2] = r[(iface2[2]-1)*3 + 2];
252 tt[3] = r[(iface2[3]-1)*3 + 2];
253 (this->*fDrawFace)(icodes, r, 4, iface2, tt);
254}
255
256////////////////////////////////////////////////////////////////////////////////
257/// Draw front surfaces of surrounding box & axes
258///
259/// \param[in] ang angle between X and Y axis
260
261void TPainter3dAlgorithms::FrontBox(Double_t ang)
262{
263 static Int_t iface1[4] = { 1, 2, 6, 5 };
264 static Int_t iface2[4] = { 2, 3, 7, 6 };
265
266 TView *view = gPad ? gPad->GetView() : nullptr;
267 if (!view) {
268 Error("FrontBox", "no TView in current pad");
269 return;
270 }
271
272 // Get corners of surrounding box
273 Double_t r[3*8], av[3*8], x[4], y[4];
274 Int_t ix1, ix2, iy1, iy2, iz1, iz2;
275 Double_t cosa = TMath::Cos(kRad*ang);
276 Double_t sina = TMath::Sin(kRad*ang);
277 view->AxisVertex(ang, av, ix1, ix2, iy1, iy2, iz1, iz2);
278 for (Int_t i = 0; i < 8; ++i) {
279 r[i*3 + 0] = av[i*3 + 0] + av[i*3 + 1]*cosa;
280 r[i*3 + 1] = av[i*3 + 1]*sina;
281 r[i*3 + 2] = av[i*3 + 2];
282 view->WCtoNDC(&r[i*3],&r[i*3]);
283 }
284
285 // Draw frame
286 SetLineColor(1);
287 SetLineStyle(1);
288 SetLineWidth(1);
289 TAttLine::Modify();
290 for (Int_t i = 0; i < 4; ++i) {
291 Int_t k = iface1[i] - 1;
292 x[i] = r[k*3 + 0];
293 y[i] = r[k*3 + 1];
294 }
295 gPad->PaintPolyLine(4, x, y);
296 for (Int_t i = 0; i < 4; ++i) {
297 Int_t k = iface2[i] - 1;
298 x[i] = r[k*3 + 0];
299 y[i] = r[k*3 + 1];
300 }
301 gPad->PaintPolyLine(4, x, y);
302}
303
304////////////////////////////////////////////////////////////////////////////////
305/// Clear screen
306
307void TPainter3dAlgorithms::ClearRaster()
308{
309 Int_t nw = (fNxrast*fNyrast + 29) / 30;
310 for (Int_t i = 0; i < nw; ++i) fRaster[i] = 0;
311 fIfrast = 0;
312}
313
314////////////////////////////////////////////////////////////////////////////////
315/// Set correspondence between function and color levels
316///
317/// \param[in] nl number of levels
318/// \param[in] fl function levels
319/// \param[in] icl colors for levels
320///
321/// \param[out] irep return code (0 OK, -1 error).
322
323void TPainter3dAlgorithms::ColorFunction(Int_t nl, Double_t *fl, Int_t *icl, Int_t &irep)
324{
325 static const char *where = "ColorFunction";
326
327 irep = 0;
328 if (nl == 0) {
329 fNlevel = 0;
330 return;
331 }
332
333 // Check parameters
334 if (nl < 0 || nl > 256) {
335 Error(where, "illegal number of levels (%d)", nl);
336 irep = -1;
337 return;
338 }
339
340 for (Int_t i = 1; i < nl; ++i) {
341 if (fl[i] <= fl[i - 1]) {
342 // Error(where, "function levels must be in increasing order");
343 irep = -1;
344 return;
345 }
346 }
347
348 for (Int_t i = 0; i < nl; ++i) {
349 if (icl[i] < 0) {
350 // Error(where, "negative color index (%d)", icl[i]);
351 irep = -1;
352 return;
353 }
354 }
355
356 // Set levels
357 fNlevel = nl;
358 for (Int_t i = 0; i < fNlevel; ++i) fFunLevel[i] = Hparam.factor*fl[i];
359 for (Int_t i = 0; i < fNlevel+1; ++i) fColorLevel[i] = icl[i];
360}
361
362////////////////////////////////////////////////////////////////////////////////
363/// Define the grid levels drawn in the background of surface and lego plots.
364/// The grid levels are aligned on the Z axis' main tick marks.
365
366void TPainter3dAlgorithms::DefineGridLevels(Int_t ndivz)
367{
368 TView *view = gPad ? gPad->GetView() : nullptr;
369 if (!view) {
370 Error("GridLevels", "no TView in current pad");
371 return;
372 }
373
374 // Find the main tick marks positions
375 Int_t nbins = 0;
376 Double_t binLow = 0, binHigh = 0, binWidth = 0;
377 Double_t *rmin = view->GetRmin();
378 Double_t *rmax = view->GetRmax();
379 if (!rmin || !rmax) return;
380 if (ndivz > 0) {
381 THLimitsFinder::Optimize(rmin[2], rmax[2], ndivz,
382 binLow, binHigh, nbins, binWidth, " ");
383 } else {
384 nbins = TMath::Abs(ndivz);
385 binLow = rmin[2];
386 binHigh = rmax[2];
387 binWidth = (binHigh - binLow)/nbins;
388 }
389
390 // Define the grid levels
391 fNlevel = nbins + 1;
392 for (Int_t i = 0; i < fNlevel; ++i) {
393 fFunLevel[i] = binLow + i*binWidth;
394 }
395}
396
397////////////////////////////////////////////////////////////////////////////////
398/// Draw face - 1st variant (2 colors: 1st for external surface, 2nd for internal)
399///
400/// \param[in] icodes set of codes for the line (not used in this method)
401/// \param[in] xyz coordinates of nodes
402/// \param[in] np number of nodes in face
403/// \param[in] iface face
404/// \param[in] t additional function defined on this face (not used in this method)
405
406void TPainter3dAlgorithms::DrawFaceMode1(Int_t *, Double_t *xyz, Int_t np, Int_t *iface, Double_t *)
407{
408 TView *view = gPad ? gPad->GetView() : nullptr;
409 if (!view) return;
410
411 // Transfer to normalised coordinates
412 Bool_t ifneg = false;
413 Double_t x[12+1] = {0}, y[12+1] = {0}, p3[3];
414 for (Int_t i = 0; i < np; ++i) {
415 Int_t k = iface[i];
416 if (k < 0) { k = -k; ifneg = true; }
417 view->WCtoNDC(&xyz[(k-1)*3], p3);
418 x[i] = p3[0]; y[i] = p3[1];
419 }
420 x[np] = x[0]; y[np] = y[0];
421
422 // Find normal
423 Double_t z = 0;
424 for (Int_t i = 0; i < np; ++i) {
425 z += y[i]*x[i+1] - x[i]*y[i+1];
426 }
427
428 // Draw face
429 SetFillColor((z > 0) ? kF3FillColor1 : kF3FillColor2);
430 SetFillStyle(1001);
431 TAttFill::Modify();
432 gPad->PaintFillArea(np, x, y);
433
434 // Draw border
435 SetLineColor(kF3LineColor);
436 TAttLine::Modify();
437 if (ifneg) {
438 for (Int_t i = 0; i < np; ++i) { // draw visible edges, skip invisible
439 if (iface[i] > 0) gPad->PaintPolyLine(2, &x[i], &y[i]);
440 }
441 } else {
442 gPad->PaintPolyLine(np+1, x, y); // all edges are visible
443 }
444}
445
446////////////////////////////////////////////////////////////////////////////////
447/// Draw face - 2nd option (fill in correspondence with function levels)
448///
449/// \param[in] icodes set of codes for the line (not used in this method)
450/// \param[in] xyz coordinates of nodes
451/// \param[in] np number of nodes
452/// \param[in] iface face
453/// \param[in] t additional function defined on this face
454
455void TPainter3dAlgorithms::DrawFaceMode2(Int_t *, Double_t *xyz, Int_t np, Int_t *iface, Double_t *t)
456{
457 TView *view = gPad ? gPad->GetView() : nullptr;
458 if (!view) return;
459
460 // Transfer to normalised coordinates
461 Double_t x[12+1] = {0}, y[12+1] = {0}, p3[3*12];
462 for (Int_t i = 0; i < np; ++i) {
463 Int_t k = iface[i];
464 view->WCtoNDC(&xyz[(k-1)*3], &p3[i*3]);
465 x[i] = p3[i*3+0]; y[i] = p3[i*3+1];
466 }
467 x[np] = x[0]; y[np] = y[0];
468
469 // Draw face
470 SetLineColor(fEdgeColor[fEdgeIdx]);
471 SetLineStyle(fEdgeStyle[fEdgeIdx]);
472 SetLineWidth(fEdgeWidth[fEdgeIdx]);
473 TAttLine::Modify();
474 if (np == 4) {
475 Double_t ttt[5] = { t[0], t[1], t[2], t[3], t[0] };
476 for (Int_t i = 0; i<3; ++i) { p3[3*4+i] = p3[i]; }
477 Int_t k1 = 0, k2 = 2;
478 Double_t z1 = (x[k1+1] - x[k1+0])*(y[k1+2] - y[k1+1]) - (y[k1+1] - y[k1+0])*(x[k1+2] - x[k1+1]);
479 Double_t z2 = (x[k2+1] - x[k2+0])*(y[k2+2] - y[k2+1]) - (y[k2+1] - y[k2+0])*(x[k2+2] - x[k2+1]);
480 if (z1 > z2) { k1 = 2; k2 = 0; }
481 FillPolygon(3, &p3[3*k1], &ttt[k1]);
482 if (fMesh == 1) { // Draw border
483 gPad->PaintPolyLine(3, &x[k1], &y[k1]);
484 }
485 FillPolygon(3, &p3[3*k2], &ttt[k2]);
486 if (fMesh == 1) { // Draw border
487 gPad->PaintPolyLine(3, &x[k2], &y[k2]);
488 if (z1*z2 <= 0) { // Draw middle line
489 x[1] = x[2]; y[1] = y[2];
490 gPad->PaintPolyLine(2, &x[0], &y[0]);
491 }
492 }
493 } else {
494 FillPolygon(np, p3, t);
495 if (fMesh == 1) { // Draw border
496 gPad->PaintPolyLine(np+1, x, y);
497 }
498 }
499}
500
501////////////////////////////////////////////////////////////////////////////////
502/// Draw face - 3rd option (draw face for stacked lego plot)
503///
504/// \param[in] icodes set of codes for the line
505/// \param[in] xyz coordinates of nodes
506/// \param[in] np number of nodes
507/// \param[in] iface face
508/// \param[in] t additional function defined on this face (not used in this method)
509
510void TPainter3dAlgorithms::DrawFaceMode3(Int_t *icodes, Double_t *xyz, Int_t np, Int_t *iface, Double_t *)
511{
512 TView *view = gPad ? gPad->GetView() : nullptr;
513 if (!view) return;
514
515 // Transfer to normalised coordinates
516 Double_t x[4+1] = {0}, y[4+1] = {0}, p3[3];
517 for (Int_t i = 0; i < np; ++i) {
518 Int_t k = iface[i];
519 view->WCtoNDC(&xyz[(k-1)*3], p3);
520 x[i] = p3[0]; y[i] = p3[1];
521 }
522 x[np] = x[0]; y[np] = y[0];
523
524 // Draw face
525 Int_t icol = 0;
526 if (icodes[3] == 6) icol = fColorTop;
527 if (icodes[3] == 5) icol = fColorBottom;
528 if (icodes[3] == 1) icol = fColorMain[icodes[2] - 1];
529 if (icodes[3] == 2) icol = fColorDark[icodes[2] - 1];
530 if (icodes[3] == 3) icol = fColorMain[icodes[2] - 1];
531 if (icodes[3] == 4) icol = fColorDark[icodes[2] - 1];
532 SetFillStyle(1001);
533 SetFillColor(icol);
534 TAttFill::Modify();
535 gPad->PaintFillArea(np, x, y);
536
537 // Draw border
538 if (fMesh) {
539 SetLineColor(fEdgeColor[fEdgeIdx]);
540 SetLineStyle(fEdgeStyle[fEdgeIdx]);
541 SetLineWidth(fEdgeWidth[fEdgeIdx]);
542 TAttLine::Modify();
543 gPad->PaintPolyLine(np+1, x, y);
544 }
545}
546
547////////////////////////////////////////////////////////////////////////////////
548/// Draw face - 1st variant for "MOVING SCREEN" algorithm (draw face with level lines)
549///
550/// \param[in] icodes set of codes for the line
551/// \param[in] xyz coordinates of nodes
552/// \param[in] np number of nodes
553/// \param[in] iface face
554/// \param[in] tt additional function defined on this face
555
556void TPainter3dAlgorithms::DrawFaceMove1(Int_t *icodes, Double_t *xyz, Int_t np,
557 Int_t *iface, Double_t *tt)
558{
559 TView *view = gPad ? gPad->GetView() : nullptr;
560 if (!view) return;
561
562 // Copy points to array
563 Double_t p3[3*12] = {0};
564 for (Int_t i = 0; i < np; ++i) {
565 Int_t k = iface[i];
566 p3[i*3 + 0] = xyz[(k-1)*3 + 0];
567 p3[i*3 + 1] = xyz[(k-1)*3 + 1];
568 p3[i*3 + 2] = xyz[(k-1)*3 + 2];
569 }
570
571 // Find level lines
572 FindLevelLines(np, p3, tt);
573
574 // Draw level lines
575 Double_t p1[3], p2[3], x[2], y[2];
576 SetLineStyle(3);
577 if (icodes[2] == 0) { // front & back boxes
578 SetLineColor(1);
579 SetLineWidth(1);
580 } else {
581 SetLineColor(fEdgeColor[fEdgeIdx]);
582 SetLineWidth(fEdgeWidth[fEdgeIdx]);
583 }
584 TAttLine::Modify();
585 for (Int_t il = 0; il < fNlines; ++il) {
586 FindVisibleDraw(&fPlines[6*il + 0], &fPlines[6*il + 3]);
587 view->WCtoNDC(&fPlines[6*il + 0], p1);
588 view->WCtoNDC(&fPlines[6*il + 3], p2);
589 Double_t xdel = p2[0] - p1[0];
590 Double_t ydel = p2[1] - p1[1];
591 for (Int_t it = 0; it < fNT; ++it) {
592 x[0] = p1[0] + xdel*fT[2*it + 0];
593 y[0] = p1[1] + ydel*fT[2*it + 0];
594 x[1] = p1[0] + xdel*fT[2*it + 1];
595 y[1] = p1[1] + ydel*fT[2*it + 1];
596 gPad->PaintPolyLine(2, x, y);
597 }
598 }
599
600 // Draw face
601 if (icodes[2] == 0) { // front & back boxes
602 SetLineColor(1);
603 SetLineStyle(1);
604 SetLineWidth(1);
605 } else {
606 SetLineColor(fEdgeColor[fEdgeIdx]);
607 SetLineStyle(fEdgeStyle[fEdgeIdx]);
608 SetLineWidth(fEdgeWidth[fEdgeIdx]);
609 }
610 TAttLine::Modify();
611 for (Int_t i = 0; i < np; ++i) {
612 Int_t i1 = i;
613 Int_t i2 = (i == np-1) ? 0 : i + 1;
614 FindVisibleDraw(&p3[i1*3], &p3[i2*3]);
615 view->WCtoNDC(&p3[i1*3], p1);
616 view->WCtoNDC(&p3[i2*3], p2);
617 Double_t xdel = p2[0] - p1[0];
618 Double_t ydel = p2[1] - p1[1];
619 for (Int_t it = 0; it < fNT; ++it) {
620 x[0] = p1[0] + xdel*fT[2*it + 0];
621 y[0] = p1[1] + ydel*fT[2*it + 0];
622 x[1] = p1[0] + xdel*fT[2*it + 1];
623 y[1] = p1[1] + ydel*fT[2*it + 1];
624 gPad->PaintPolyLine(2, x, y);
625 }
626 }
627
628 // Modify screen
629 for (Int_t i = 0; i < np; ++i) {
630 Int_t i1 = i;
631 Int_t i2 = (i == np-1) ? 0 : i + 1;
632 ModifyScreen(&p3[i1*3], &p3[i2*3]);
633 }
634}
635
636////////////////////////////////////////////////////////////////////////////////
637/// Draw face - 2nd variant for "MOVING SCREEN" algorithm (draw face for stacked lego plot)
638///
639/// \param[in] icodes set of codes for the line
640/// \param[in] xyz coordinates of nodes
641/// \param[in] np number of nodes
642/// \param[in] iface face
643/// \param[in] tt additional function defined on this face (not used in this method)
644
645void TPainter3dAlgorithms::DrawFaceMove2(Int_t *icodes, Double_t *xyz, Int_t np, Int_t *iface, Double_t *)
646{
647 TView *view = gPad ? gPad->GetView() : nullptr;
648 if (!view) return;
649
650 // Copy points to array
651 Double_t p3[3*12];
652 for (Int_t i = 0; i < np; ++i) {
653 Int_t k = iface[i];
654 p3[i*3 + 0] = xyz[(k-1)*3 + 0];
655 p3[i*3 + 1] = xyz[(k-1)*3 + 1];
656 p3[i*3 + 2] = xyz[(k-1)*3 + 2];
657 }
658
659 // Draw face
660 Double_t p1[3], p2[3], x[2], y[2];
661 if (icodes[2] == 0) { // front & back boxes
662 SetLineColor(1);
663 SetLineStyle(1);
664 SetLineWidth(1);
665 } else {
666 SetLineColor(fEdgeColor[fEdgeIdx]);
667 SetLineStyle(fEdgeStyle[fEdgeIdx]);
668 SetLineWidth(fEdgeWidth[fEdgeIdx]);
669 }
670 TAttLine::Modify();
671 for (Int_t i = 0; i < np; ++i) {
672 Int_t i1 = i;
673 Int_t i2 = (i == np-1) ? 0 : i + 1;
674 FindVisibleDraw(&p3[i1*3], &p3[i2*3]);
675 view->WCtoNDC(&p3[i1*3], p1);
676 view->WCtoNDC(&p3[i2*3], p2);
677 Double_t xdel = p2[0] - p1[0];
678 Double_t ydel = p2[1] - p1[1];
679 for (Int_t it = 0; it < fNT; ++it) {
680 x[0] = p1[0] + xdel*fT[2*it + 0];
681 y[0] = p1[1] + ydel*fT[2*it + 0];
682 x[1] = p1[0] + xdel*fT[2*it + 1];
683 y[1] = p1[1] + ydel*fT[2*it + 1];
684 gPad->PaintPolyLine(2, x, y);
685 }
686 }
687
688 // Modify screen
689 for (Int_t i = 0; i < np; ++i) {
690 Int_t i1 = i;
691 Int_t i2 = (i == np-1) ? 0 : i + 1;
692 ModifyScreen(&p3[i1*3], &p3[i2*3]);
693 }
694}
695
696////////////////////////////////////////////////////////////////////////////////
697/// Draw face - 3rd variant for "MOVING SCREEN" algorithm (draw level lines only)
698///
699/// \param[in] icodes set of codes for the line
700/// \param[in] xyz coordinates of nodes
701/// \param[in] np number of nodes
702/// \param[in] iface face
703/// \param[in] tt additional function defined on this face
704
705void TPainter3dAlgorithms::DrawFaceMove3(Int_t *icodes, Double_t *xyz, Int_t np,
706 Int_t *iface, Double_t *tt)
707{
708 TView *view = gPad ? gPad->GetView() : nullptr;
709 if (!view) return;
710
711 // Set graphics attributes
712 if (icodes[2] == 0) { // frame
713 SetLineColor(1);
714 SetLineStyle(1);
715 SetLineWidth(1);
716 } else {
717 SetLineColor(fEdgeColor[fEdgeIdx]);
718 SetLineStyle(fEdgeStyle[fEdgeIdx]);
719 SetLineWidth(fEdgeWidth[fEdgeIdx]);
720 }
721 TAttLine::Modify();
722
723 // Copy points to array
724 Double_t p3[3*12] = {0}, ttt[12] = {0};
725 for (Int_t i = 0; i < np; ++i) {
726 Int_t k = iface[i];
727 p3[i*3 + 0] = xyz[(k-1)*3 + 0];
728 p3[i*3 + 1] = xyz[(k-1)*3 + 1];
729 p3[i*3 + 2] = xyz[(k-1)*3 + 2];
730 ttt[i] = tt[i];
731 }
732
733 // Subdivide quadrilateral in two triangles
734 Int_t npol[2] = { np, 0 }; // number of vertices in sub-polygons
735 Int_t ipol[2] = { 0, 0 }; // first vertices in sub-polygons
736 if (np == 4 && icodes[2] != 0) {
737 p3[4*3 + 0] = p3[0];
738 p3[4*3 + 1] = p3[1];
739 p3[4*3 + 2] = p3[2];
740 ttt[4] = tt[0];
741 npol[0] = 3; npol[1] = 3;
742 ipol[0] = 0; ipol[1] = 2;
743 }
744
745 Double_t p1[3], p2[3], x[2], y[2];
746 for (Int_t kpol = 0; kpol < 2; ++kpol) {
747 if (npol[kpol] == 0) continue;
748 Int_t nv = npol[kpol];
749 Int_t iv = ipol[kpol];
750
751 // Find level lines
752 FindLevelLines(nv, &p3[3*iv], &ttt[iv]);
753
754 // Draw level lines
755 for (Int_t il = 0; il < fNlines; ++il) {
756 FindVisibleDraw(&fPlines[6*il + 0], &fPlines[6*il + 3]);
757 view->WCtoNDC(&fPlines[6*il + 0], p1);
758 view->WCtoNDC(&fPlines[6*il + 3], p2);
759 Double_t xdel = p2[0] - p1[0];
760 Double_t ydel = p2[1] - p1[1];
761 for (Int_t it = 0; it < fNT; ++it) {
762 x[0] = p1[0] + xdel*fT[2*it + 0];
763 y[0] = p1[1] + ydel*fT[2*it + 0];
764 x[1] = p1[0] + xdel*fT[2*it + 1];
765 y[1] = p1[1] + ydel*fT[2*it + 1];
766 gPad->PaintPolyLine(2, x, y);
767 }
768 }
769 }
770
771 // Modify screen
772 for (Int_t i = 0; i < np; ++i) {
773 Int_t i1 = i;
774 Int_t i2 = (i == np - 1) ? 0 : i1 + 1;
775 ModifyScreen(&p3[i1*3], &p3[i2*3]);
776 }
777}
778
779////////////////////////////////////////////////////////////////////////////////
780/// Draw level lines without hidden line removal
781///
782/// \param[in] icodes set of codes for the line
783/// \param[in] xyz coordinates of nodes
784/// \param[in] np number of nodes
785/// \param[in] iface face
786/// \param[in] tt additional function defined on this face
787
788void TPainter3dAlgorithms::DrawLevelLines(Int_t *icodes, Double_t *xyz, Int_t np,
789 Int_t *iface, Double_t *tt)
790{
791 TView *view = gPad ? gPad->GetView() : nullptr;
792 if (!view) return;
793
794 // Set graphics attributes
795 if (icodes[2] == 0) { // frame
796 SetLineColor(1);
797 SetLineStyle(1);
798 SetLineWidth(1);
799 } else {
800 SetLineColor(fEdgeColor[fEdgeIdx]);
801 SetLineStyle(fEdgeStyle[fEdgeIdx]);
802 SetLineWidth(fEdgeWidth[fEdgeIdx]);
803 }
804 TAttLine::Modify();
805
806 // Copy points to array
807 Double_t p3[3*12] = {0}, ttt[12] = {0};
808 for (Int_t i = 0; i < np; ++i) {
809 Int_t k = iface[i];
810 p3[i*3 + 0] = xyz[(k-1)*3 + 0];
811 p3[i*3 + 1] = xyz[(k-1)*3 + 1];
812 p3[i*3 + 2] = xyz[(k-1)*3 + 2];
813 ttt[i] = tt[i];
814 }
815
816 // Subdivide quadrilateral in two triangles
817 Int_t npol[2] = { np, 0 }; // number of vertices in sub-polygons
818 Int_t ipol[2] = { 0, 0 }; // first vertices in sub-polygons
819 if (np == 4 && icodes[2] != 0) {
820 p3[4*3 + 0] = p3[0];
821 p3[4*3 + 1] = p3[1];
822 p3[4*3 + 2] = p3[2];
823 ttt[4] = tt[0];
824 npol[0] = 3; npol[1] = 3;
825 ipol[0] = 0; ipol[1] = 2;
826 }
827
828 Double_t p1[3], p2[3], x[2], y[2];
829 for (Int_t kpol = 0; kpol < 2; ++kpol) {
830 if (npol[kpol] == 0) continue;
831 Int_t nv = npol[kpol];
832 Int_t iv = ipol[kpol];
833
834 // Find level lines
835 FindLevelLines(nv, &p3[3*iv], &ttt[iv]);
836
837 // Draw level lines
838 for (Int_t il = 0; il < fNlines; ++il) {
839 view->WCtoNDC(&fPlines[6*il + 0], p1);
840 view->WCtoNDC(&fPlines[6*il + 3], p2);
841 x[0] = p1[0]; y[0] = p1[1];
842 x[1] = p2[0]; y[1] = p2[1];
843 gPad->PaintPolyLine(2, x, y);
844 }
845 }
846}
847
848////////////////////////////////////////////////////////////////////////////////
849/// Draw face - 1st variant for "RASTER SCREEN" algorithm (draw face with level lines)
850///
851/// \param[in] icodes set of codes for the line
852/// \param[in] xyz coordinates of nodes
853/// \param[in] np number of nodes
854/// \param[in] iface face
855/// \param[in] tt additional function defined on this face
856
857void TPainter3dAlgorithms::DrawFaceRaster1(Int_t *icodes, Double_t *xyz, Int_t np, Int_t *iface, Double_t *tt)
858{
859 TView *view = gPad ? gPad->GetView() : nullptr;
860 if (!view) return;
861
862 // Copy vertices to array
863 Double_t p3[3*12] = {0}, pp[2*12] = {0};
864 for (Int_t i = 0; i < np; ++i) {
865 Int_t k = iface[i];
866 if (k < 0) k = -k;
867 p3[i*3 + 0] = xyz[(k-1)*3 + 0];
868 p3[i*3 + 1] = xyz[(k-1)*3 + 1];
869 p3[i*3 + 2] = xyz[(k-1)*3 + 2];
870 Double_t p[3];
871 view->WCtoNDC(&p3[i*3], p);
872 pp[2*i + 0] = p[0];
873 pp[2*i + 1] = p[1];
874 }
875
876 // Find level lines
877 FindLevelLines(np, p3, tt);
878
879 // Draw level lines
880 Double_t p1[3], p2[3], x[2], y[2];
881 SetLineStyle(3);
882 if (icodes[2] == 0) { // front & back boxes
883 SetLineColor(1);
884 SetLineWidth(1);
885 } else {
886 SetLineColor(fEdgeColor[fEdgeIdx]);
887 SetLineWidth(fEdgeWidth[fEdgeIdx]);
888 }
889 TAttLine::Modify();
890 for (Int_t il = 0; il < fNlines; ++il) {
891 view->WCtoNDC(&fPlines[6*il + 0], p1);
892 view->WCtoNDC(&fPlines[6*il + 3], p2);
893 FindVisibleLine(p1, p2, 100, fNT, fT);
894 Double_t xdel = p2[0] - p1[0];
895 Double_t ydel = p2[1] - p1[1];
896 for (Int_t it = 0; it < fNT; ++it) {
897 x[0] = p1[0] + xdel*fT[2*it + 0];
898 y[0] = p1[1] + ydel*fT[2*it + 0];
899 x[1] = p1[0] + xdel*fT[2*it + 1];
900 y[1] = p1[1] + ydel*fT[2*it + 1];
901 gPad->PaintPolyLine(2, x, y);
902 }
903 }
904
905 // Draw face
906 if (icodes[2] == 0) { // front & back boxes
907 SetLineColor(1);
908 SetLineStyle(1);
909 SetLineWidth(1);
910 } else {
911 SetLineColor(fEdgeColor[fEdgeIdx]);
912 SetLineStyle(fEdgeStyle[fEdgeIdx]);
913 SetLineWidth(fEdgeWidth[fEdgeIdx]);
914 }
915 TAttLine::Modify();
916 for (Int_t i = 0; i < np; ++i) {
917 if (iface[i] < 0) continue;
918 Int_t i1 = i;
919 Int_t i2 = (i == np-1) ? 0 : i + 1;
920 FindVisibleLine(&pp[2*i1], &pp[2*i2], 100, fNT, fT);
921 Double_t xdel = pp[2*i2 + 0] - pp[2*i1 + 0];
922 Double_t ydel = pp[2*i2 + 1] - pp[2*i1 + 1];
923 for (Int_t it = 0; it < fNT; ++it) {
924 x[0] = pp[2*i1 + 0] + xdel*fT[2*it + 0];
925 y[0] = pp[2*i1 + 1] + ydel*fT[2*it + 0];
926 x[1] = pp[2*i1 + 0] + xdel*fT[2*it + 1];
927 y[1] = pp[2*i1 + 1] + ydel*fT[2*it + 1];
928 gPad->PaintPolyLine(2, x, y);
929 }
930 }
931
932 // Modify raster screen
933 FillPolygonBorder(np, pp);
934}
935
936////////////////////////////////////////////////////////////////////////////////
937/// Draw face - 2nd variant for "RASTER SCREEN" algorithm (draw face for stacked lego plot)
938///
939/// \param[in] icodes set of codes for the line (not used in this method)
940/// \param[in] xyz coordinates of nodes
941/// \param[in] np number of nodes
942/// \param[in] iface face
943/// \param[in] tt additional function defined on this face (not used in this method)
944
945void TPainter3dAlgorithms::DrawFaceRaster2(Int_t *, Double_t *xyz, Int_t np, Int_t *iface, Double_t *)
946{
947 TView *view = gPad ? gPad->GetView() : nullptr;
948 if (!view) return;
949
950 // Copy vertices to array
951 Double_t x[2], y[2], pp[2*12];
952 for (Int_t i = 0; i < np; ++i) {
953 Int_t k = iface[i];
954 if (k < 0) k = -k;
955 Double_t p[3];
956 view->WCtoNDC(&xyz[(k-1)*3], p);
957 pp[2*i + 0] = p[0];
958 pp[2*i + 1] = p[1];
959 }
960
961 // Draw face
962 SetLineColor(fEdgeColor[fEdgeIdx]);
963 SetLineStyle(fEdgeStyle[fEdgeIdx]);
964 SetLineWidth(fEdgeWidth[fEdgeIdx]);
965 TAttLine::Modify();
966 for (Int_t i = 0; i < np; ++i) {
967 if (iface[i] < 0) continue;
968 Int_t i1 = i;
969 Int_t i2 = (i == np-1) ? 0 : i + 1;
970 FindVisibleLine(&pp[2*i1], &pp[2*i2], 100, fNT, fT);
971 Double_t xdel = pp[2*i2 + 0] - pp[2*i1 + 0];
972 Double_t ydel = pp[2*i2 + 1] - pp[2*i1 + 1];
973 for (Int_t it = 0; it < fNT; ++it) {
974 x[0] = pp[2*i1 + 0] + xdel*fT[2*it + 0];
975 y[0] = pp[2*i1 + 1] + ydel*fT[2*it + 0];
976 x[1] = pp[2*i1 + 0] + xdel*fT[2*it + 1];
977 y[1] = pp[2*i1 + 1] + ydel*fT[2*it + 1];
978 gPad->PaintPolyLine(2, x, y);
979 }
980 }
981
982 // Modify raster screen
983 FillPolygonBorder(np, pp);
984}
985
986////////////////////////////////////////////////////////////////////////////////
987/// Fill polygon with function values at vertexes
988///
989/// \param[in] n number of vertexes
990/// \param[in] p polygon
991/// \param[in] f function values at nodes
992///
993/// Errors:
994/// - illegal number of vertexes in polygon
995/// - illegal call of FillPolygon: no levels
996
997void TPainter3dAlgorithms::FillPolygon(Int_t n, Double_t *p, Double_t *f)
998{
999 Int_t ilev, i, k, icol, i1, i2, nl, np;
1000 Double_t fmin, fmax;
1001 Double_t x[12], y[12], f1, f2;
1002 Double_t p3[36] /* was [3][12] */;
1003 Double_t funmin, funmax;
1004
1005 /* Parameter adjustments */
1006 --f;
1007 p -= 4;
1008
1009 if (n < 3) {
1010 Error("FillPolygon", "illegal number of vertices in polygon (%d)", n);
1011 return;
1012 }
1013
1014 if (fNlevel == 0) {
1015 // Illegal call of FillPolygon: no levels
1016 return;
1017 }
1018 np = n;
1019 nl = fNlevel;
1020 if (nl < 0) nl = -nl;
1021 fmin = f[1];
1022 fmax = f[1];
1023 for (i = 2; i <= np; ++i) {
1024 if (fmin > f[i]) fmin = f[i];
1025 if (fmax < f[i]) fmax = f[i];
1026 }
1027 funmin = fFunLevel[0] - 1;
1028 if (fmin < funmin) funmin = fmin - 1;
1029 funmax = fFunLevel[nl - 1] + 1;
1030 if (fmax > funmax) funmax = fmax + 1;
1031
1032 // F I N D A N D D R A W S U B P O L Y G O N S
1033 f2 = funmin;
1034 for (ilev = 1; ilev <= nl+1; ++ilev) {
1035 // S E T L E V E L L I M I T S
1036 f1 = f2;
1037 if (ilev == nl + 1) f2 = funmax;
1038 else f2 = fFunLevel[ilev - 1];
1039 if (fmax < f1) return;
1040 if (fmin > f2) continue;
1041 // F I N D S U B P O L Y G O N
1042 k = 0;
1043 for (i = 1; i <= np; ++i) {
1044 i1 = i;
1045 i2 = i + 1;
1046 if (i == np) i2 = 1;
1047 FindPartEdge(&p[i1*3 + 1], &p[i2*3 + 1], f[i1], f[i2], f1, f2, k, p3);
1048 }
1049 // D R A W S U B P O L Y G O N
1050 if (k < 3) continue;
1051 for (i = 1; i <= k; ++i) {
1052 x[i-1] = p3[i*3-3];
1053 y[i-1] = p3[i*3-2];
1054 if (TMath::IsNaN(x[i-1]) || TMath::IsNaN(y[i-1])) return;
1055 }
1056 if (ilev==1) {
1057 icol=gPad->GetFillColor();
1058 } else {
1059 icol = fColorLevel[ilev - 2];
1060 }
1061 SetFillColor(icol);
1062 SetFillStyle(1001);
1063 TAttFill::Modify();
1064 gPad->PaintFillArea(k, x, y);
1065 }
1066}
1067
1068////////////////////////////////////////////////////////////////////////////////
1069/// Fill a polygon including border ("RASTER SCREEN")
1070///
1071/// \param[in] nn number of polygon nodes
1072/// \param[in] xy polygon nodes
1073
1074void TPainter3dAlgorithms::FillPolygonBorder(Int_t nn, Double_t *xy)
1075{
1076 Int_t kbit, nbit, step, ymin, ymax, test[kLmax], xcur[kLmax], xnex[kLmax],
1077 i, j, k, n, ibase, t, x, y, xscan[24] /* was [2][kLmax] */,
1078 yscan, x1[kLmax+2], y1[kLmax+2], x2[kLmax+2], y2[kLmax+2],
1079 ib, nb, dx, dy, iw, nx, xx, yy, signdx, nstart, xx1, xx2, nxa, nxb;
1080
1081 // T R A N S F E R T O S C R E E N C O O R D I N A T E S
1082 /* Parameter adjustments */
1083 xy -= 3;
1084
1085 if (fIfrast) return;
1086
1087 n = nn;
1088 x1[0] = 0;
1089 y1[0] = 0;
1090 for (i = 1; i <= n; ++i) {
1091 x1[i - 1] = Int_t(fNxrast*((xy[2*i + 1] - fXrast) /fDXrast) - 0.01);
1092 y1[i - 1] = Int_t(fNyrast*((xy[2*i + 2] - fYrast) /fDYrast) - 0.01);
1093 }
1094 x1[n] = x1[0];
1095 y1[n] = y1[0];
1096
1097 // F I N D Y - M I N A N D Y - M A X
1098 // S E T R I G H T E D G E O R I E N T A T I O N
1099 ymin = y1[0];
1100 ymax = y1[0];
1101 for (i = 1; i <= n; ++i) {
1102 if (ymin > y1[i - 1]) ymin = y1[i - 1];
1103 if (ymax < y1[i - 1]) ymax = y1[i - 1];
1104 if (y1[i - 1] <= y1[i]) {x2[i - 1] = x1[i]; y2[i - 1] = y1[i];}
1105 else {
1106 x2[i - 1] = x1[i - 1];
1107 y2[i - 1] = y1[i - 1];
1108 x1[i - 1] = x1[i];
1109 y1[i - 1] = y1[i];
1110 }
1111 }
1112 if (ymin >= fNyrast) return;
1113 if (ymax < 0) return;
1114 if (ymax >= fNyrast) ymax = fNyrast - 1;
1115
1116 // S O R T L I N E S
1117 for (i = 1; i < n; ++i) {
1118 if (y1[i] >= y1[i - 1]) continue;
1119 y = y1[i];
1120 k = 1;
1121 for (j = i - 1; j >= 1; --j) {
1122 if (y < y1[j - 1]) continue;
1123 k = j + 1;
1124 break;
1125 }
1126 x = x1[i];
1127 xx = x2[i];
1128 yy = y2[i];
1129 for (j = i; j >= k; --j) {
1130 x1[j] = x1[j - 1];
1131 y1[j] = y1[j - 1];
1132 x2[j] = x2[j - 1];
1133 y2[j] = y2[j - 1];
1134 }
1135 x1[k - 1] = x;
1136 y1[k - 1] = y;
1137 x2[k - 1] = xx;
1138 y2[k - 1] = yy;
1139 }
1140
1141 // S E T I N I T I A L V A L U E S
1142 for (i = 1; i <= n; ++i) {
1143 xcur[i - 1] = x1[i - 1];
1144 dy = y2[i - 1] - y1[i - 1];
1145 dx = x2[i - 1] - x1[i - 1];
1146 signdx = 1;
1147 if (dx < 0) signdx = -1;
1148 if (dx < 0) dx = -dx;
1149 if (dx <= dy) {
1150 t = -(dy + 1) / 2 + dx;
1151 if (t < 0) {
1152 test[i - 1] = t;
1153 xnex[i - 1] = xcur[i - 1];
1154 } else {
1155 test[i - 1] = t - dy;
1156 xnex[i - 1] = xcur[i - 1] + signdx;
1157 }
1158 } else if (dy != 0) {
1159 step = (dx - 1) / (dy + dy) + 1;
1160 test[i - 1] = step*dy - (dx + 1) / 2 - dx;
1161 xnex[i - 1] = xcur[i - 1] + signdx*step;
1162 }
1163 }
1164
1165 // L O O P O N S C A N L I N E S
1166 nstart = 1;
1167 for (yscan = ymin; yscan <= ymax; ++yscan) {
1168 nx = 0;
1169 nxa = 0;
1170 nxb = kLmax + 1;
1171 for (i = nstart; i <= n; ++i) {
1172 if (y1[i - 1] > yscan) goto L500;
1173 if (y2[i - 1] <= yscan) {
1174 if (i == nstart) ++nstart;
1175 if (y2[i - 1] != yscan)continue;
1176 --nxb;
1177 if (x2[i - 1] >= xcur[i - 1]) {
1178 xscan[2*nxb - 2] = xcur[i - 1];
1179 xscan[2*nxb - 1] = x2[i - 1];
1180 } else {
1181 xscan[2*nxb - 2] = x2[i - 1];
1182 xscan[2*nxb - 1] = xcur[i - 1];
1183 }
1184 continue;
1185 }
1186
1187 // S T O R E C U R R E N T X
1188 // P R E P A R E X F O R N E X T S C A N - L I N E
1189 ++nxa;
1190 dy = y2[i - 1] - y1[i - 1];
1191 dx = x2[i - 1] - x1[i - 1];
1192 if (dx >= 0) {
1193 signdx = 1;
1194 xscan[2*nxa - 2] = xcur[i - 1];
1195 xscan[2*nxa - 1] = xnex[i - 1];
1196 if (xscan[2*nxa - 2] != xscan[2*nxa - 1]) {
1197 --xscan[2*nxa - 1];
1198 }
1199 } else {
1200 dx = -dx;
1201 signdx = -1;
1202 xscan[2*nxa - 2] = xnex[i - 1];
1203 xscan[2*nxa - 1] = xcur[i - 1];
1204 if (xscan[2*nxa - 2] != xscan[2*nxa - 1]) {
1205 ++xscan[2*nxa - 2];
1206 }
1207 }
1208 xcur[i - 1] = xnex[i - 1];
1209 if (dx <= dy) {
1210 test[i - 1] += dx;
1211 if (test[i - 1] < 0) continue;
1212 test[i - 1] -= dy;
1213 xnex[i - 1] += signdx;
1214 continue;
1215 }
1216 step = dx / dy;
1217 t = test[i - 1] + step*dy;
1218 if (t >= 0) {
1219 test[i - 1] = t - dx;
1220 xnex[i - 1] += signdx*step;
1221 } else {
1222 test[i - 1] = t + dy - dx;
1223 xnex[i - 1] += signdx*(step + 1);
1224 }
1225 }
1226
1227 // S O R T P O I N T S A L O N G X
1228L500:
1229 if (yscan < 0) continue;
1230 ibase = yscan*fNxrast;
1231 if (nxa >= 2) {
1232 for (i = 1; i < nxa; ++i) {
1233 for (j = i; j >= 1; --j) {
1234 if (xscan[2*j] >= xscan[2*j - 2]) continue;
1235 x = xscan[2*j];
1236 xscan[2*j] = xscan[2*j - 2];
1237 xscan[2*j - 2] = x;
1238 x = xscan[2*j - 1];
1239 xscan[2*j + 1] = xscan[2*j - 1];
1240 xscan[2*j - 1] = x;
1241 }
1242 }
1243 for (i = 1; i <= nxa; i += 2) {
1244 ++nx;
1245 xscan[2*nx - 2] = xscan[2*i - 2];
1246 x = xscan[2*i + 1];
1247 if (xscan[2*i - 1] > x) x = xscan[2*i - 1];
1248 xscan[2*nx - 1] = x;
1249 }
1250 }
1251 if (nxb <= kLmax) {
1252 for (i = nxb; i <= kLmax; ++i) {
1253 ++nx;
1254 xscan[2*nx - 2] = xscan[2*i - 2];
1255 xscan[2*nx - 1] = xscan[2*i - 1];
1256 }
1257 }
1258 // C O N C A T E N A T E A N D F I L L
1259 while (nx) {
1260 xx1 = xscan[2*nx - 2];
1261 xx2 = xscan[2*nx - 1];
1262 --nx;
1263 k = 1;
1264 while (k <= nx) {
1265 if ((xscan[2*k - 2] <= xx2 + 1) && (xscan[2*k - 1] >= xx1 - 1)) {
1266 if (xscan[2*k - 2] < xx1) xx1 = xscan[2*k - 2];
1267 if (xscan[2*k - 1] > xx2) xx2 = xscan[2*k - 1];
1268 xscan[2*k - 2] = xscan[2*nx - 2];
1269 xscan[2*k - 1] = xscan[2*nx - 1];
1270 --nx;
1271 } else ++k;
1272 }
1273 if (xx1 < 0) xx1 = 0;
1274 if (xx2 >= fNxrast) xx2 = fNxrast - 1;
1275 nbit = xx2 - xx1 + 1;
1276 kbit = ibase + xx1;
1277 iw = kbit / 30;
1278 ib = kbit - iw*30 + 1;
1279 iw = iw + 1;
1280 nb = 30 - ib + 1;
1281 if (nb > nbit) nb = nbit;
1282 fRaster[iw - 1] = fRaster[iw - 1] | fMask[fJmask[nb - 1] + ib - 1];
1283 nbit -= nb;
1284 if (nbit) {
1285 while(nbit > 30) {
1286 fRaster[iw] = fMask[464];
1287 ++iw;
1288 nbit += -30;
1289 }
1290 fRaster[iw] = fRaster[iw] | fMask[fJmask[nbit - 1]];
1291 ++iw;
1292 }
1293 }
1294 }
1295}
1296
1297////////////////////////////////////////////////////////////////////////////////
1298/// Find level lines for face
1299///
1300/// \param[in] np number of nodes
1301/// \param[in] f face
1302/// \param[in] t additional function
1303///
1304/// Error: number of points for line not equal 2
1305
1306void TPainter3dAlgorithms::FindLevelLines(Int_t np, Double_t *f, Double_t *t)
1307{
1308 fNlines = 0;
1309 if (fNlevel == 0) return;
1310 Int_t nl = TMath::Abs(fNlevel);
1311
1312 // Find Tmin and Tmax
1313 Double_t tmin = t[0];
1314 Double_t tmax = t[0];
1315 for (Int_t i = 1; i < np; ++i) {
1316 if (t[i] < tmin) tmin = t[i];
1317 if (t[i] > tmax) tmax = t[i];
1318 }
1319 if (tmin >= fFunLevel[nl - 1]) return;
1320 if (tmax <= fFunLevel[0]) return;
1321
1322 // Find level lines
1323 for (Int_t il = 1; il <= nl; ++il) {
1324 if (tmin >= fFunLevel[il - 1]) continue;
1325 if (tmax < fFunLevel[il - 1]) return;
1326 if (fNlines >= 200) return;
1327 fNlines++;
1328 fLevelLine[fNlines - 1] = il;
1329 Int_t kp = 0;
1330 for (Int_t i = 0; i < np; ++i) {
1331 Int_t i1 = i;
1332 Int_t i2 = (i == np-1) ? 0 : i+1;
1333 Double_t d1 = t[i1] - fFunLevel[il - 1];
1334 Double_t d2 = t[i2] - fFunLevel[il - 1];
1335 if (d1 == 0) d1 = 1e-99;
1336 if (d2 == 0) d2 = 1e-99;
1337 if (d1*d2 > 0) continue;
1338
1339 // find point
1340 kp++;
1341 d1 /= t[i2] - t[i1];
1342 d2 /= t[i2] - t[i1];
1343 fPlines[(kp + 2*fNlines)*3 - 9] = d2*f[i1*3 + 0] - d1*f[i2*3 + 0];
1344 fPlines[(kp + 2*fNlines)*3 - 8] = d2*f[i1*3 + 1] - d1*f[i2*3 + 1];
1345 fPlines[(kp + 2*fNlines)*3 - 7] = d2*f[i1*3 + 2] - d1*f[i2*3 + 2];
1346 if (kp == 2) break;
1347 }
1348 if (kp != 2) {
1349 Error("FindLevelLines", "number of points for line not equal 2");
1350 fNlines--;
1351 }
1352 }
1353}
1354
1355////////////////////////////////////////////////////////////////////////////////
1356/// Find part of edge where function defined on this edge has value from
1357/// `fmin` to `fmax`
1358///
1359/// \param[in] p1 1st point
1360/// \param[in] p2 2nd point
1361/// \param[in] f1 function value at 1st point
1362/// \param[in] f2 function value at 2nd point
1363/// \param[in] fmin min value of layer
1364/// \param[in] fmax max value of layer
1365///
1366/// \param[out] kpp current number of point
1367/// \param[out] pp coordinates of new face
1368
1369void TPainter3dAlgorithms::FindPartEdge(Double_t *p1, Double_t *p2, Double_t f1,
1370 Double_t f2, Double_t fmin,
1371 Double_t fmax, Int_t &kpp, Double_t *pp)
1372{
1373 Double_t d1, d2;
1374 Int_t k1, k2, kk;
1375
1376 /* Parameter adjustments */
1377 pp -= 4;
1378 --p2;
1379 --p1;
1380
1381 k1 = 0;
1382 if (f1 < fmin) k1 = -2;
1383 if (f1 == fmin) k1 = -1;
1384 if (f1 == fmax) k1 = 1;
1385 if (f1 > fmax) k1 = 2;
1386 k2 = 0;
1387 if (f2 < fmin) k2 = -2;
1388 if (f2 == fmin) k2 = -1;
1389 if (f2 == fmax) k2 = 1;
1390 if (f2 > fmax) k2 = 2;
1391 kk = (k1 + 2)*5 + (k2 + 2) + 1;
1392
1393 // K2: -2 -1 0 +1 +2
1394 // K1: -2 -1 0 +1 +2
1395 switch ((int)kk) {
1396 case 1: return;
1397 case 2: return;
1398 case 3: goto L200;
1399 case 4: goto L200;
1400 case 5: goto L600;
1401 case 6: goto L100;
1402 case 7: goto L100;
1403 case 8: goto L100;
1404 case 9: goto L100;
1405 case 10: goto L500;
1406 case 11: goto L400;
1407 case 12: goto L100;
1408 case 13: goto L100;
1409 case 14: goto L100;
1410 case 15: goto L500;
1411 case 16: goto L400;
1412 case 17: goto L100;
1413 case 18: goto L100;
1414 case 19: goto L100;
1415 case 20: goto L100;
1416 case 21: goto L700;
1417 case 22: goto L300;
1418 case 23: goto L300;
1419 case 24: return;
1420 case 25: return;
1421 }
1422
1423 // 1 - S T P O I N T
1424L100:
1425 ++kpp;
1426 pp[kpp*3 + 1] = p1[1];
1427 pp[kpp*3 + 2] = p1[2];
1428 pp[kpp*3 + 3] = p1[3];
1429 return;
1430
1431 // I N T E R S E C T I O N W I T H Fmin
1432L200:
1433 ++kpp;
1434 d1 = (fmin - f1) / (f1 - f2);
1435 d2 = (fmin - f2) / (f1 - f2);
1436 pp[kpp*3 + 1] = d2*p1[1] - d1*p2[1];
1437 pp[kpp*3 + 2] = d2*p1[2] - d1*p2[2];
1438 pp[kpp*3 + 3] = d2*p1[3] - d1*p2[3];
1439 return;
1440
1441 // I N T E R S E C T I O N W I T H Fmax
1442L300:
1443 ++kpp;
1444 d1 = (fmax - f1) / (f1 - f2);
1445 d2 = (fmax - f2) / (f1 - f2);
1446 pp[kpp*3 + 1] = d2*p1[1] - d1*p2[1];
1447 pp[kpp*3 + 2] = d2*p1[2] - d1*p2[2];
1448 pp[kpp*3 + 3] = d2*p1[3] - d1*p2[3];
1449 return;
1450
1451 // 1 - S T P O I N T, I N T E R S E C T I O N WITH Fmin
1452L400:
1453 ++kpp;
1454 pp[kpp*3 + 1] = p1[1];
1455 pp[kpp*3 + 2] = p1[2];
1456 pp[kpp*3 + 3] = p1[3];
1457 ++kpp;
1458 d1 = (fmin - f1) / (f1 - f2);
1459 d2 = (fmin - f2) / (f1 - f2);
1460 pp[kpp*3 + 1] = d2*p1[1] - d1*p2[1];
1461 pp[kpp*3 + 2] = d2*p1[2] - d1*p2[2];
1462 pp[kpp*3 + 3] = d2*p1[3] - d1*p2[3];
1463 return;
1464
1465 // 1 - S T P O I N T, I N T E R S E C T I O N WITH Fmax
1466L500:
1467 ++kpp;
1468 pp[kpp*3 + 1] = p1[1];
1469 pp[kpp*3 + 2] = p1[2];
1470 pp[kpp*3 + 3] = p1[3];
1471 ++kpp;
1472 d1 = (fmax - f1) / (f1 - f2);
1473 d2 = (fmax - f2) / (f1 - f2);
1474 pp[kpp*3 + 1] = d2*p1[1] - d1*p2[1];
1475 pp[kpp*3 + 2] = d2*p1[2] - d1*p2[2];
1476 pp[kpp*3 + 3] = d2*p1[3] - d1*p2[3];
1477 return;
1478
1479 // I N T E R S E C T I O N W I T H Fmin, Fmax
1480L600:
1481 ++kpp;
1482 d1 = (fmin - f1) / (f1 - f2);
1483 d2 = (fmin - f2) / (f1 - f2);
1484 pp[kpp*3 + 1] = d2*p1[1] - d1*p2[1];
1485 pp[kpp*3 + 2] = d2*p1[2] - d1*p2[2];
1486 pp[kpp*3 + 3] = d2*p1[3] - d1*p2[3];
1487 ++kpp;
1488 d1 = (fmax - f1) / (f1 - f2);
1489 d2 = (fmax - f2) / (f1 - f2);
1490 pp[kpp*3 + 1] = d2*p1[1] - d1*p2[1];
1491 pp[kpp*3 + 2] = d2*p1[2] - d1*p2[2];
1492 pp[kpp*3 + 3] = d2*p1[3] - d1*p2[3];
1493 return;
1494
1495 // I N T E R S E C T I O N W I T H Fmax, Fmin
1496L700:
1497 ++kpp;
1498 d1 = (fmax - f1) / (f1 - f2);
1499 d2 = (fmax - f2) / (f1 - f2);
1500 pp[kpp*3 + 1] = d2*p1[1] - d1*p2[1];
1501 pp[kpp*3 + 2] = d2*p1[2] - d1*p2[2];
1502 pp[kpp*3 + 3] = d2*p1[3] - d1*p2[3];
1503 ++kpp;
1504 d1 = (fmin - f1) / (f1 - f2);
1505 d2 = (fmin - f2) / (f1 - f2);
1506 pp[kpp*3 + 1] = d2*p1[1] - d1*p2[1];
1507 pp[kpp*3 + 2] = d2*p1[2] - d1*p2[2];
1508 pp[kpp*3 + 3] = d2*p1[3] - d1*p2[3];
1509}
1510
1511////////////////////////////////////////////////////////////////////////////////
1512/// Find visible parts of line (draw line)
1513///
1514/// \param[in] r1 1-st point of the line
1515/// \param[in] r2 2-nd point of the line
1516
1517void TPainter3dAlgorithms::FindVisibleDraw(Double_t *r1, Double_t *r2)
1518{
1519 Double_t yy1u, yy2u;
1520 Int_t i, icase, i1, i2, icase1, icase2, iv, ifback;
1521 Double_t x1, x2, y1, y2, z1, z2, dd, di;
1522 Double_t dt, dy;
1523 Double_t tt, uu, ww, yy, yy1, yy2, yy1d, yy2d;
1524 Double_t *tn = nullptr;
1525 const Double_t kEpsil = 1.e-6;
1526 /* Parameter adjustments */
1527 --r2;
1528 --r1;
1529 TView *view = gPad ? gPad->GetView() : nullptr;
1530
1531 if (view) {
1532 tn = view->GetTN();
1533 if (tn) {
1534 x1 = tn[0]*r1[1] + tn[1]*r1[2] + tn[2]*r1[3] + tn[3];
1535 x2 = tn[0]*r2[1] + tn[1]*r2[2] + tn[2]*r2[3] + tn[3];
1536 y1 = tn[4]*r1[1] + tn[5]*r1[2] + tn[6]*r1[3] + tn[7];
1537 y2 = tn[4]*r2[1] + tn[5]*r2[2] + tn[6]*r2[3] + tn[7];
1538 z1 = tn[8]*r1[1] + tn[9]*r1[2] + tn[10]*r1[3] + tn[11];
1539 z2 = tn[8]*r2[1] + tn[9]*r2[2] + tn[10]*r2[3] + tn[11];
1540 } else {
1541 Error("FindVisibleDraw", "invalid TView in current pad");
1542 return;
1543 }
1544 } else {
1545 Error("FindVisibleDraw", "no TView in current pad");
1546 return;
1547 }
1548
1549 ifback = 0;
1550 if (x1 >= x2) {
1551 ifback = 1;
1552 ww = x1;
1553 x1 = x2;
1554 x2 = ww;
1555 ww = y1;
1556 y1 = y2;
1557 y2 = ww;
1558 ww = z1;
1559 z1 = z2;
1560 z2 = ww;
1561 }
1562 fNT = 0;
1563 i1 = Int_t((x1 - fX0) / fDX) + 15;
1564 i2 = Int_t((x2 - fX0) / fDX) + 15;
1565 x1 = fX0 + (i1 - 1)*fDX;
1566 x2 = fX0 + (i2 - 1)*fDX;
1567 if (i1 != i2) {
1568
1569 // F I N D V I S I B L E P A R T S O F T H E L I N E
1570 di = (Double_t) (i2 - i1);
1571 dy = (y2 - y1) / di;
1572 dt = 1 / di;
1573 iv = -1;
1574 for (i = i1; i <= i2 - 1; ++i) {
1575 yy1 = y1 + dy*(i - i1);
1576 yy2 = yy1 + dy;
1577 yy1u = yy1 - fU[2*i - 2];
1578 yy1d = yy1 - fD[2*i - 2];
1579 yy2u = yy2 - fU[2*i - 1];
1580 yy2d = yy2 - fD[2*i - 1];
1581 tt = dt*(i - i1);
1582 // A N A L I Z E L E F T S I D E
1583 icase1 = 1;
1584 if (yy1u > kEpsil) icase1 = 0;
1585 if (yy1d < -kEpsil) icase1 = 2;
1586 if ((icase1 == 0 || icase1 == 2) && iv <= 0) {
1587 iv = 1;
1588 ++fNT;
1589 fT[2*fNT - 2] = tt;
1590 }
1591 if (icase1 == 1 && iv >= 0) {
1592 iv = -1;
1593 fT[2*fNT - 1] = tt;
1594 }
1595 // A N A L I Z E R I G H T S I D E
1596 icase2 = 1;
1597 if (yy2u > kEpsil) icase2 = 0;
1598 if (yy2d < -kEpsil) icase2 = 2;
1599 icase = icase1*3 + icase2;
1600 if (icase == 1) {
1601 iv = -1;
1602 fT[2*fNT - 1] = tt + dt*(yy1u / (yy1u - yy2u));
1603 }
1604 if (icase == 2) {
1605 fT[2*fNT - 1] = tt + dt*(yy1u / (yy1u - yy2u));
1606 ++fNT;
1607 fT[2*fNT - 2] = tt + dt*(yy1d / (yy1d - yy2d));
1608 }
1609 if (icase == 3) {
1610 iv = 1;
1611 ++fNT;
1612 fT[2*fNT - 2] = tt + dt*(yy1u / (yy1u - yy2u));
1613 }
1614 if (icase == 5) {
1615 iv = 1;
1616 ++fNT;
1617 fT[2*fNT - 2] = tt + dt*(yy1d / (yy1d - yy2d));
1618 }
1619 if (icase == 6) {
1620 fT[2*fNT - 1] = tt + dt*(yy1d / (yy1d - yy2d));
1621 ++fNT;
1622 fT[2*fNT - 2] = tt + dt*(yy1u / (yy1u - yy2u));
1623 }
1624 if (icase == 7) {
1625 iv = -1;
1626 fT[2*fNT - 1] = tt + dt*(yy1d / (yy1d - yy2d));
1627 }
1628 if (fNT + 1 >= 100) break;
1629 }
1630 if (iv > 0) fT[2*fNT - 1] = 1;
1631 } else {
1632
1633 // V E R T I C A L L I N E
1634 fNT = 1;
1635 fT[0] = 0;
1636 fT[1] = 1;
1637 if (y2 <= y1) {
1638 if (y2 == y1) { fNT = 0; return;}
1639 ifback = 1 - ifback;
1640 yy = y1;
1641 y1 = y2;
1642 y2 = yy;
1643 }
1644 uu = fU[2*i1 - 2];
1645 dd = fD[2*i1 - 2];
1646 if (i1 != 1) {
1647 if (uu < fU[2*i1 - 3]) uu = fU[2*i1 - 3];
1648 if (dd > fD[2*i1 - 3]) dd = fD[2*i1 - 3];
1649 }
1650 // F I N D V I S I B L E P A R T O F L I N E
1651 if (y1 < uu && y2 > dd) {
1652 if (y1 >= dd && y2 <= uu) {fNT = 0; return;}
1653 fNT = 0;
1654 if (dd > y1) {
1655 ++fNT;
1656 fT[2*fNT - 2] = 0;
1657 fT[2*fNT - 1] = (dd - y1) / (y2 - y1);
1658 }
1659 if (uu < y2) {
1660 ++fNT;
1661 fT[2*fNT - 2] = (uu - y1) / (y2 - y1);
1662 fT[2*fNT - 1] = 1;
1663 }
1664 }
1665 }
1666
1667 if (ifback == 0) return;
1668 if (fNT == 0) return;
1669 for (i = 1; i <= fNT; ++i) {
1670 fT[2*i - 2] = 1 - fT[2*i - 2];
1671 fT[2*i - 1] = 1 - fT[2*i - 1];
1672 }
1673}
1674
1675////////////////////////////////////////////////////////////////////////////////
1676/// Find visible part of a line ("RASTER SCREEN")
1677///
1678/// \param[in] p1 1st point of the line
1679/// \param[in] p2 2nd point of the line
1680/// \param[in] ntmax max allowed number of visible segments
1681///
1682/// \param[out] nt number of visible segments of the line
1683/// \param[out] t visible segments
1684
1685void TPainter3dAlgorithms::FindVisibleLine(Double_t *p1, Double_t *p2, Int_t ntmax, Int_t &nt, Double_t *t)
1686{
1687 Double_t ddtt;
1688 Double_t tcur;
1689 Int_t i, incrx, ivis, x1, y1, x2, y2, ib, kb, dx, dy, iw, ix, iy, ifinve, dx2, dy2;
1690 Double_t t1, t2;
1691 Double_t dt;
1692 Double_t tt;
1693 /* Parameter adjustments */
1694 t -= 3;
1695 --p2;
1696 --p1;
1697
1698 if (fIfrast) {
1699 nt = 1;
1700 t[3] = 0;
1701 t[4] = 1;
1702 return;
1703 }
1704 x1 = Int_t(fNxrast*((p1[1] - fXrast) / fDXrast) - 0.01);
1705 y1 = Int_t(fNyrast*((p1[2] - fYrast) / fDYrast) - 0.01);
1706 x2 = Int_t(fNxrast*((p2[1] - fXrast) / fDXrast) - 0.01);
1707 y2 = Int_t(fNyrast*((p2[2] - fYrast) / fDYrast) - 0.01);
1708 ifinve = 0;
1709 if (y1 > y2) {
1710 ifinve = 1;
1711 iw = x1;
1712 x1 = x2;
1713 x2 = iw;
1714 iw = y1;
1715 y1 = y2;
1716 y2 = iw;
1717 }
1718 nt = 0;
1719 ivis = 0;
1720 if (y1 >= fNyrast) return;
1721 if (y2 < 0) return;
1722 if (x1 >= fNxrast && x2 >= fNxrast) return;
1723 if (x1 < 0 && x2 < 0) return;
1724
1725 // S E T I N I T I A L V A L U E S
1726 incrx = 1;
1727 dx = x2 - x1;
1728 if (dx < 0) {
1729 dx = -dx;
1730 incrx = -1;
1731 }
1732 dy = y2 - y1;
1733 dx2 = dx + dx;
1734 dy2 = dy + dy;
1735 if (dy > dx) goto L200;
1736
1737 // D X . G T . D Y
1738 dt = 1./ (Double_t)(dx + 1.);
1739 ddtt = dt*(float).5;
1740 tcur = -(Double_t)dt;
1741 tt = (Double_t) (-(dx + dy2));
1742 iy = y1;
1743 kb = iy*fNxrast + x1 - incrx;
1744 for (ix = x1; incrx < 0 ? ix >= x2 : ix <= x2; ix += incrx) {
1745 kb += incrx;
1746 tcur += dt;
1747 tt += dy2;
1748 if (tt >= 0) {
1749 ++iy;
1750 tt -= dx2;
1751 kb += fNxrast;
1752 }
1753 if (iy < 0) goto L110;
1754 if (iy >= fNyrast) goto L110;
1755 if (ix < 0) goto L110;
1756 if (ix >= fNxrast) goto L110;
1757 iw = kb / 30;
1758 ib = kb - iw*30 + 1;
1759 if (fRaster[iw] & fMask[ib - 1]) goto L110;
1760 if (ivis > 0) continue;
1761 ivis = 1;
1762 ++nt;
1763 t[2*nt + 1] = tcur;
1764 continue;
1765L110:
1766 if (ivis == 0) continue;
1767 ivis = 0;
1768 t[2*nt + 2] = tcur;
1769 if (nt == ntmax) goto L300;
1770 }
1771 if (ivis > 0) t[2*nt + 2] = tcur + dt + ddtt;
1772 goto L300;
1773
1774 // D Y . G T . D X
1775L200:
1776 dt = 1. / (Double_t)(dy + 1.);
1777 ddtt = dt*(float).5;
1778 tcur = -(Double_t)dt;
1779 tt = (Double_t) (-(dy + dx2));
1780 ix = x1;
1781 if (y2 >= fNyrast) y2 = fNyrast - 1;
1782 kb = (y1 - 1)*fNxrast + ix;
1783 for (iy = y1; iy <= y2; ++iy) {
1784 kb += fNxrast;
1785 tcur += dt;
1786 tt += dx2;
1787 if (tt >= 0) {
1788 ix += incrx;
1789 tt -= dy2;
1790 kb += incrx;
1791 }
1792 if (iy < 0) goto L210;
1793 if (ix < 0) goto L210;
1794 if (ix >= fNxrast) goto L210;
1795 iw = kb / 30;
1796 ib = kb - iw*30 + 1;
1797 if (fRaster[iw] & fMask[ib - 1]) goto L210;
1798 if (ivis > 0) continue;
1799 ivis = 1;
1800 ++nt;
1801 t[2*nt + 1] = tcur;
1802 continue;
1803L210:
1804 if (ivis == 0) continue;
1805 ivis = 0;
1806 t[2*nt + 2] = tcur;
1807 if (nt == ntmax) goto L300;
1808 }
1809 if (ivis > 0) t[2*nt + 2] = tcur + dt;
1810
1811 // C H E C K D I R E C T I O N O F P A R A M E T E R
1812L300:
1813 if (nt == 0) return;
1814 dt *= 1.1;
1815 if (t[3] <= dt) t[3] = 0;
1816 if (t[2*nt + 2] >= 1 - dt) t[2*nt + 2] = 1;
1817 if (ifinve == 0) return;
1818 for (i = 1; i <= nt; ++i) {
1819 t1 = t[2*i + 1];
1820 t2 = t[2*i + 2];
1821 t[2*i + 1] = 1 - t2;
1822 t[2*i + 2] = 1 - t1;
1823 }
1824}
1825
1826////////////////////////////////////////////////////////////////////////////////
1827/// Find part of surface with luminosity in the corners. This method is used for
1828/// Gouraud shading
1829
1830void TPainter3dAlgorithms::GouraudFunction(Int_t ia, Int_t ib, Double_t *face, Double_t *t)
1831{
1832 Int_t iphi;
1833 static Double_t f[108]; /* was [3][4][3][3] */
1834 Int_t i, j, k;
1835 Double_t r, s, x[36]; /* was [4][3][3] */
1836 Double_t y[36]; /* was [4][3][3] */
1837 Double_t z[36]; /* was [4][3][3] */
1838 Int_t incrx[3], incry[3];
1839
1840 Double_t x1, x2, y1, y2, z1, z2, th, an[27]; /* was [3][3][3] */
1841 Double_t bn[12]; /* was [3][2][2] */
1842
1843 Double_t rad;
1844 Double_t phi;
1845 Int_t ixt, iyt;
1846
1847 /* Parameter adjustments */
1848 --t;
1849 face -= 4;
1850
1851 iphi = 1;
1852 rad = TMath::ATan(1) * (float)4 / (float)180;
1853
1854 // Find real cell indexes
1855 ixt = ia + Hparam.xfirst - 1;
1856 iyt = ib + Hparam.yfirst - 1;
1857
1858 // Find increments of neighboring cells
1859 incrx[0] = -1;
1860 incrx[1] = 0;
1861 incrx[2] = 1;
1862 if (ixt == 1) incrx[0] = 0;
1863 if (ixt == Hparam.xlast - 1) incrx[2] = 0;
1864 incry[0] = -1;
1865 incry[1] = 0;
1866 incry[2] = 1;
1867 if (iyt == 1) incry[0] = 0;
1868 if (iyt == Hparam.ylast - 1) incry[2] = 0;
1869
1870 // Find neighboring faces
1871 Int_t i1, i2;
1872 for (j = 1; j <= 3; ++j) {
1873 for (i = 1; i <= 3; ++i) {
1874 i1 = ia + incrx[i - 1];
1875 i2 = ib + incry[j - 1];
1876 SurfaceFunction(i1, i2, &f[(((i + j*3) << 2) + 1)*3 - 51], &t[1]);
1877 }
1878 }
1879
1880 // Set face
1881 for (k = 1; k <= 4; ++k) {
1882 for (i = 1; i <= 3; ++i) {
1883 face[i + k*3] = f[i + (k + 32)*3 - 52];
1884 }
1885 }
1886
1887 // Find coordinates and normales
1888 for (j = 1; j <= 3; ++j) {
1889 for (i = 1; i <= 3; ++i) {
1890 for (k = 1; k <= 4; ++k) {
1891 if (Hoption.System == kPOLAR) {
1892 phi = f[iphi + (k + ((i + j*3) << 2))*3 - 52]*rad;
1893 r = f[3 - iphi + (k + ((i + j*3) << 2))*3 - 52];
1894 x[k + ((i + j*3) << 2) - 17] = r * TMath::Cos(phi);
1895 y[k + ((i + j*3) << 2) - 17] = r * TMath::Sin(phi);
1896 z[k + ((i + j*3) << 2) - 17] = f[(k + ((i + j*3) << 2))*3 - 49];
1897 } else if (Hoption.System == kCYLINDRICAL) {
1898 phi = f[iphi + (k + ((i + j*3) << 2))*3 - 52]*rad;
1899 r = f[(k + ((i + j*3) << 2))*3 - 49];
1900 x[k + ((i + j*3) << 2) - 17] = r*TMath::Cos(phi);
1901 y[k + ((i + j*3) << 2) - 17] = r*TMath::Sin(phi);
1902 z[k + ((i + j*3) << 2) - 17] = f[3 - iphi + (k + ((i + j*3) << 2))*3 - 52];
1903 } else if (Hoption.System == kSPHERICAL) {
1904 phi = f[iphi + (k + ((i + j*3) << 2))*3 - 52]*rad;
1905 th = f[3 - iphi + (k + ((i + j*3) << 2))*3 - 52]*rad;
1906 r = f[(k + ((i + j*3) << 2))*3 - 49];
1907 x[k + ((i + j*3) << 2) - 17] = r*TMath::Sin(th)*TMath::Cos(phi);
1908 y[k + ((i + j*3) << 2) - 17] = r*TMath::Sin(th)*TMath::Sin(phi);
1909 z[k + ((i + j*3) << 2) - 17] = r*TMath::Cos(th);
1910 } else if (Hoption.System == kRAPIDITY) {
1911 phi = f[iphi + (k + ((i + j*3) << 2))*3 - 52]*rad;
1912 th = f[3 - iphi + (k + ((i + j*3) << 2))*3 - 52]*rad;
1913 r = f[(k + ((i + j*3) << 2))*3 - 49];
1914 x[k + ((i + j*3) << 2) - 17] = r*TMath::Cos(phi);
1915 y[k + ((i + j*3) << 2) - 17] = r*TMath::Sin(phi);
1916 z[k + ((i + j*3) << 2) - 17] = r*TMath::Cos(th) / TMath::Sin(th);
1917 } else {
1918 x[k + ((i + j*3) << 2) - 17] = f[(k + ((i + j*3) << 2))*3 - 51];
1919 y[k + ((i + j*3) << 2) - 17] = f[(k + ((i + j*3) << 2))*3 - 50];
1920 z[k + ((i + j*3) << 2) - 17] = f[(k + ((i + j*3) << 2))*3 - 49];
1921 }
1922 }
1923 x1 = x[((i + j*3) << 2) - 14] - x[((i + j*3) << 2) - 16];
1924 x2 = x[((i + j*3) << 2) - 13] - x[((i + j*3) << 2) - 15];
1925 y1 = y[((i + j*3) << 2) - 14] - y[((i + j*3) << 2) - 16];
1926 y2 = y[((i + j*3) << 2) - 13] - y[((i + j*3) << 2) - 15];
1927 z1 = z[((i + j*3) << 2) - 14] - z[((i + j*3) << 2) - 16];
1928 z2 = z[((i + j*3) << 2) - 13] - z[((i + j*3) << 2) - 15];
1929 an[(i + j*3)*3 - 12] = y1*z2 - y2*z1;
1930 an[(i + j*3)*3 - 11] = z1*x2 - z2*x1;
1931 an[(i + j*3)*3 - 10] = x1*y2 - x2*y1;
1932 s = TMath::Sqrt(an[(i + j*3)*3 - 12]*an[(i + j*3)*3 - 12] + an[
1933 (i + j*3)*3 - 11]*an[(i + j*3)*3 - 11] + an[(i
1934 + j*3)*3 - 10]*an[(i + j*3)*3 - 10]);
1935
1936 an[(i + j*3)*3 - 12] /= s;
1937 an[(i + j*3)*3 - 11] /= s;
1938 an[(i + j*3)*3 - 10] /= s;
1939 }
1940 }
1941
1942 // Find average normals
1943 for (j = 1; j <= 2; ++j) {
1944 for (i = 1; i <= 2; ++i) {
1945 for (k = 1; k <= 3; ++k) {
1946 bn[k + (i + 2*j)*3 - 10] = an[k + (i + j*3)*3 - 13]
1947 + an[k + (i + 1 + j*3)*3 - 13] + an[k + (i + 1 +
1948 (j + 1)*3)*3 - 13] + an[k + (i + (j + 1)*3)*3 - 13];
1949 }
1950 }
1951 }
1952
1953 TView *view = gPad ? gPad->GetView() : nullptr;
1954
1955 // Set luminosity
1956 Luminosity(view, bn, t[1]);
1957 Luminosity(view, &bn[3], t[2]);
1958 Luminosity(view, &bn[9], t[3]);
1959 Luminosity(view, &bn[6], t[4]);
1960}
1961
1962////////////////////////////////////////////////////////////////////////////////
1963/// Initialize "MOVING SCREEN" method
1964///
1965/// \param[in] xmin left boundary
1966/// \param[in] xmax right boundary
1967
1968void TPainter3dAlgorithms::InitMoveScreen(Double_t xmin, Double_t xmax)
1969{
1970 const Double_t VERY_BIG = 9e+99;
1971 fX0 = xmin;
1972 fDX = (xmax - xmin) / NumOfSlices;
1973 for (Int_t i = 0; i < NumOfSlices; ++i) {
1974 fU[2*i + 0] = -VERY_BIG;
1975 fU[2*i + 1] = -VERY_BIG;
1976 fD[2*i + 0] = VERY_BIG;
1977 fD[2*i + 1] = VERY_BIG;
1978 }
1979}
1980
1981////////////////////////////////////////////////////////////////////////////////
1982/// Initialize hidden lines removal algorithm (RASTER SCREEN)
1983///
1984/// \param[in] xmin Xmin in the normalized coordinate system
1985/// \param[in] ymin Ymin in the normalized coordinate system
1986/// \param[in] xmax Xmax in the normalized coordinate system
1987/// \param[in] ymax Ymax in the normalized coordinate system
1988/// \param[in] nx number of pixels along X
1989/// \param[in] ny number of pixels along Y
1990
1991void TPainter3dAlgorithms::InitRaster(Double_t xmin, Double_t ymin, Double_t xmax, Double_t ymax, Int_t nx, Int_t ny )
1992{
1993 Int_t i, j, k, ib, nb;
1994
1995 fNxrast = nx;
1996 fNyrast = ny;
1997 fXrast = xmin;
1998 fDXrast = xmax - xmin;
1999 fYrast = ymin;
2000 fDYrast = ymax - ymin;
2001
2002 // Create buffer for raster
2003 Int_t buffersize = nx*ny/30 + 1;
2004 fRaster.resize(buffersize);
2005
2006 // S E T M A S K S
2007 k = 0;
2008 Int_t pow2 = 1;
2009 for (i = 1; i <= 30; ++i) {
2010 fJmask[i - 1] = k;
2011 k = k + 30 - i + 1;
2012 fMask[i - 1] = pow2;
2013 pow2 *= 2;
2014 }
2015 j = 30;
2016 for (nb = 2; nb <= 30; ++nb) {
2017 for (ib = 1; ib <= 30 - nb + 1; ++ib) {
2018 k = 0;
2019 for (i = ib; i <= ib + nb - 1; ++i) k = k | fMask[i - 1];
2020 ++j;
2021 fMask[j - 1] = k;
2022 }
2023 }
2024
2025 // C L E A R R A S T E R S C R E E N
2026 ClearRaster();
2027}
2028
2029////////////////////////////////////////////////////////////////////////////////
2030/// Service function for Legos
2031
2032void TPainter3dAlgorithms::LegoFunction(Int_t ia, Int_t ib, Int_t &nv, Double_t *ab, Double_t *vv, Double_t *t)
2033{
2034 Int_t i, j, ixt, iyt;
2035 Double_t yval1l, yval2l;
2036 Double_t xlab1l, xlab2l, ylab1l, ylab2l;
2037 Double_t rinrad = gStyle->GetLegoInnerR();
2038 Double_t dangle = 10; //Delta angle for Rapidity option
2039
2040 /* Parameter adjustments */
2041 t -= 5;
2042 --vv;
2043 ab -= 3;
2044
2045 ixt = ia + Hparam.xfirst - 1;
2046 iyt = ib + Hparam.yfirst - 1;
2047
2048 // Compute the cell position in cartesian coordinates
2049 // and compute the LOG if necessary
2050 Double_t xwid = gCurrentHist->GetXaxis()->GetBinWidth(ixt);
2051 Double_t ywid = gCurrentHist->GetYaxis()->GetBinWidth(iyt);
2052 ab[3] = gCurrentHist->GetXaxis()->GetBinLowEdge(ixt) + xwid*Hparam.baroffset;
2053 ab[4] = gCurrentHist->GetYaxis()->GetBinLowEdge(iyt) + ywid*Hparam.baroffset;
2054 ab[5] = ab[3] + xwid*Hparam.barwidth;
2055 ab[8] = ab[4] + ywid*Hparam.barwidth;
2056
2057 if (Hoption.Logx) {
2058 if (ab[3] > 0) ab[3] = TMath::Log10(ab[3]);
2059 else ab[3] = Hparam.xmin;
2060 if (ab[5] > 0) ab[5] = TMath::Log10(ab[5]);
2061 else ab[5] = Hparam.xmin;
2062 }
2063 // xval1l = Hparam.xmin;
2064 // xval2l = Hparam.xmax;
2065 if (Hoption.Logy) {
2066 if (ab[4] > 0) ab[4] = TMath::Log10(ab[4]);
2067 else ab[4] = Hparam.ymin;
2068 if (ab[8] > 0) ab[8] = TMath::Log10(ab[8]);
2069 else ab[8] = Hparam.ymin;
2070 }
2071 yval1l = Hparam.ymin;
2072 yval2l = Hparam.ymax;
2073
2074 if (ab[3] < Hparam.xmin) ab[3] = Hparam.xmin;
2075 if (ab[4] < Hparam.ymin) ab[4] = Hparam.ymin;
2076 if (ab[5] > Hparam.xmax) ab[5] = Hparam.xmax;
2077 if (ab[8] > Hparam.ymax) ab[8] = Hparam.ymax;
2078 if (ab[5] < Hparam.xmin) ab[5] = Hparam.xmin;
2079 if (ab[8] < Hparam.ymin) ab[8] = Hparam.ymin;
2080
2081 xlab1l = gCurrentHist->GetXaxis()->GetXmin();
2082 xlab2l = gCurrentHist->GetXaxis()->GetXmax();
2083 if (Hoption.Logx) {
2084 if (xlab2l>0) {
2085 if (xlab1l>0) xlab1l = TMath::Log10(xlab1l);
2086 else xlab1l = TMath::Log10(0.001*xlab2l);
2087 xlab2l = TMath::Log10(xlab2l);
2088 }
2089 }
2090 ylab1l = gCurrentHist->GetYaxis()->GetXmin();
2091 ylab2l = gCurrentHist->GetYaxis()->GetXmax();
2092 if (Hoption.Logy) {
2093 if (ylab2l>0) {
2094 if (ylab1l>0) ylab1l = TMath::Log10(ylab1l);
2095 else ylab1l = TMath::Log10(0.001*ylab2l);
2096 ylab2l = TMath::Log10(ylab2l);
2097 }
2098 }
2099
2100 // Transform the cell position in the required coordinate system
2101 if (Hoption.System == kPOLAR) {
2102 ab[3] = 360*(ab[3] - xlab1l) / (xlab2l - xlab1l);
2103 ab[5] = 360*(ab[5] - xlab1l) / (xlab2l - xlab1l);
2104 ab[4] = (ab[4] - yval1l) / (yval2l - yval1l);
2105 ab[8] = (ab[8] - yval1l) / (yval2l - yval1l);
2106 } else if (Hoption.System == kCYLINDRICAL) {
2107 ab[3] = 360*(ab[3] - xlab1l) / (xlab2l - xlab1l);
2108 ab[5] = 360*(ab[5] - xlab1l) / (xlab2l - xlab1l);
2109 } else if (Hoption.System == kSPHERICAL) {
2110 ab[3] = 360*(ab[3] - xlab1l) / (xlab2l - xlab1l);
2111 ab[5] = 360*(ab[5] - xlab1l) / (xlab2l - xlab1l);
2112 ab[4] = 180*(ab[4] - ylab1l) / (ylab2l - ylab1l);
2113 ab[8] = 180*(ab[8] - ylab1l) / (ylab2l - ylab1l);
2114 } else if (Hoption.System == kRAPIDITY) {
2115 ab[3] = 360*(ab[3] - xlab1l) / (xlab2l - xlab1l);
2116 ab[5] = 360*(ab[5] - xlab1l) / (xlab2l - xlab1l);
2117 ab[4] = (180 - dangle*2)*(ab[4] - ylab1l) / (ylab2l - ylab1l) + dangle;
2118 ab[8] = (180 - dangle*2)*(ab[8] - ylab1l) / (ylab2l - ylab1l) + dangle;
2119 }
2120
2121 // Complete the cell coordinates
2122 ab[6] = ab[4];
2123 ab[7] = ab[5];
2124 ab[9] = ab[3];
2125 ab[10] = ab[8];
2126
2127 // Get the content of the table, and loop on the
2128 // stack if necessary.
2129 vv[1] = Hparam.zmin;
2130 vv[2] = Hparam.factor*gCurrentHist->GetBinContent(ixt, iyt);
2131
2132 // In linear scale, 3D boxes all start from 0.
2133 if (Hparam.zmin<0 && !Hoption.Logz && Hoption.MinimumZero) {
2134 if (vv[2]<0) {
2135 vv[1] = vv[2];
2136 vv[2] = 0;
2137 } else {
2138 vv[1] = 0;
2139 }
2140 }
2141
2142 TList *stack = gCurrentHist->GetPainter()->GetStack();
2143 Int_t nids = 0; //not yet implemented
2144 if (stack) nids = stack->GetSize();
2145 if (nids) {
2146 for (i = 2; i <= nids + 1; ++i) {
2147 TH1 *hid = (TH1*)stack->At(i-2);
2148 vv[i + 1] = Hparam.factor*hid->GetBinContent(ixt, iyt) + vv[i];
2149 vv[i + 1] = TMath::Max(Hparam.zmin, vv[i + 1]);
2150 //vv[i + 1] = TMath::Min(Hparam.zmax, vv[i + 1]);
2151 }
2152 }
2153
2154 nv = nids + 2;
2155 for (i = 2; i <= nv; ++i) {
2156 if (Hoption.Logz) {
2157 if (vv[i] > 0)
2158 vv[i] = TMath::Max(Hparam.zmin, (Double_t)TMath::Log10(vv[i]));
2159 else
2160 vv[i] = Hparam.zmin;
2161 vv[i] = TMath::Min(vv[i], Hparam.zmax);
2162 } else {
2163 vv[i] = TMath::Max(Hparam.zmin, vv[i]);
2164 vv[i] = TMath::Min(Hparam.zmax, vv[i]);
2165 }
2166 }
2167
2168 if (!Hoption.Logz) {
2169 i = 3;
2170 while (i <= nv) {
2171 if (vv[i] < vv[i - 1]) {
2172 vv[i - 1] = vv[i];
2173 i = 3;
2174 continue;
2175 }
2176 ++i;
2177 }
2178 }
2179
2180 // For cylindrical, spherical and pseudo-rapidity, the content
2181 // is mapped onto the radius
2182 if (Hoption.System == kCYLINDRICAL || Hoption.System == kSPHERICAL || Hoption.System == kRAPIDITY) {
2183 for (i = 1; i <= nv; ++i) {
2184 vv[i] = (1 - rinrad)*((vv[i] - Hparam.zmin) /
2185 (Hparam.zmax - Hparam.zmin)) + rinrad;
2186 }
2187 }
2188
2189 for (i = 1; i <= nv; ++i) {
2190 for (j = 1; j <= 4; ++j) t[j + (i << 2)] = vv[i];
2191 }
2192}
2193
2194////////////////////////////////////////////////////////////////////////////////
2195/// Draw stack of lego-plots in cartesian coordinates
2196///
2197/// \param[in] ang angle between X ang Y (not used in this method)
2198/// \param[in] nx number of cells along X
2199/// \param[in] ny number of cells along Y
2200/// \param[in] chopt specific options
2201///
2202/// - `chopt` = 'BF' from BACK to FRONT
2203/// - `chopt` = 'FB' from FRONT to BACK
2204
2205void TPainter3dAlgorithms::LegoCartesian(Double_t, Int_t nx, Int_t ny, const char *chopt)
2206{
2207 Int_t icodes[4], iface[4];
2208 Double_t xy[4*2], xyz[8*3], tface[4];
2209 Int_t firstStackNumberDrawn=-1 ; // necessary to compute fColorBottom when the 0 option is set and when the stack is seen from below (bottomview, theta<0.)
2210
2211 TView *view = gPad ? gPad->GetView() : nullptr;
2212 if (!view) {
2213 Error("LegoCartesian", "no TView in current pad");
2214 return;
2215 }
2216 Double_t *tnorm = view->GetTnorm();
2217 if (!tnorm) return;
2218
2219 // Allocate v and tt arrays
2220 Int_t vSize = fNStack+2;
2221 std::vector<Double_t> v(vSize), tt(4*vSize);
2222
2223 // Define order of drawing
2224 Int_t incrx = (tnorm[8] < 0.) ? -1 : +1;
2225 Int_t incry = (tnorm[9] < 0.) ? -1 : +1;
2226 if (*chopt != 'B' && *chopt != 'b') { // front to back
2227 incrx = -incrx; incry = -incry;
2228 }
2229 Int_t ix1 = (incrx == +1) ? 1 : nx;
2230 Int_t iy1 = (incry == +1) ? 1 : ny;
2231 Int_t ix2 = (incrx == +1) ? nx : 1;
2232 Int_t iy2 = (incry == +1) ? ny : 1;
2233
2234 // Find visibility of sides
2235 Double_t zn;
2236 Int_t ivis[6] = { 0,0,0,0,0,0 };
2237 view->FindNormal(0, 1, 0, zn);
2238 if (zn < 0) ivis[0] = 1;
2239 if (zn > 0) ivis[2] = 1;
2240 view->FindNormal(1, 0, 0, zn);
2241 if (zn > 0) ivis[1] = 1;
2242 if (zn < 0) ivis[3] = 1;
2243 view->FindNormal(0, 0, 1, zn);
2244 if (zn > 0) ivis[5] = 1;
2245 if (zn < 0) ivis[4] = 1;
2246
2247 // Draw stack of lego-plots
2248 Int_t nv = 0;
2249 THistPainter *painter = (THistPainter*)gCurrentHist->GetPainter();
2250 for (Int_t iy = iy1; iy != iy2+incry; iy += incry) {
2251 for (Int_t ix = ix1; ix != ix2+incrx; ix += incrx) {
2252 if (!painter->IsInside(ix,iy)) continue;
2253 (this->*fLegoFunction)(ix, iy, nv, xy, v.data(), tt.data());
2254 if (nv < 2 || nv > vSize) continue;
2255 if (Hoption.Zero) {
2256 Double_t total_content = 0;
2257 for (Int_t iv = 1; iv < nv; ++iv) { total_content += v[iv]; }
2258 if (total_content <= Hparam.zmin) continue;
2259 }
2260 icodes[0] = ix;
2261 icodes[1] = iy;
2262 for (Int_t i = 1; i <= 4; ++i) {
2263 xyz[i*3 - 3] = xy[2*i - 2];
2264 xyz[i*3 - 2] = xy[2*i - 1];
2265 xyz[(i + 4)*3 - 3] = xyz[i*3 - 3];
2266 xyz[(i + 4)*3 - 2] = xyz[i*3 - 2];
2267 }
2268 // Draw stack
2269 firstStackNumberDrawn = -1;
2270 for (Int_t iv = 1; iv < nv; ++iv) {
2271 for (Int_t i = 1; i <= 4; ++i) {
2272 xyz[i*3 - 1] = v[iv - 1];
2273 xyz[(i + 4)*3 - 1] = v[iv];
2274 }
2275 if (v[iv - 1] == v[iv]) continue;
2276 icodes[2] = iv;
2277 for (Int_t i = 1; i <= 4; ++i) {
2278 if (ivis[i - 1] == 0) continue;
2279 Int_t k1 = i;
2280 Int_t k2 = i + 1;
2281 if (i == 4) k2 = 1;
2282 icodes[3] = k1;
2283 iface[0] = k1;
2284 iface[1] = k2;
2285 iface[2] = k2 + 4;
2286 iface[3] = k1 + 4;
2287 tface[0] = tt[k1 + (iv << 2) - 5];
2288 tface[1] = tt[k2 + (iv << 2) - 5];
2289 tface[2] = tt[k2 + ((iv + 1) << 2) - 5];
2290 tface[3] = tt[k1 + ((iv + 1) << 2) - 5];
2291 fEdgeIdx = iv-1;
2292 (this->*fDrawFace)(icodes, xyz, 4, iface, tface);
2293 }
2294 if ( firstStackNumberDrawn==-1 ) firstStackNumberDrawn = fEdgeIdx;
2295 }
2296 // Draw bottom face
2297 if (ivis[4] > 0) {
2298 icodes[2] = 1;
2299 icodes[3] = 5;
2300 for (Int_t i = 1; i <= 4; ++i) {
2301 xyz[i*3 - 1] = v[0];
2302 iface[i - 1] = 5 - i;
2303 tface[i - 1] = tt[5 - i - 1];
2304 }
2305 if (!Hoption.Zero) fEdgeIdx = 0;
2306 else {
2307 fEdgeIdx = firstStackNumberDrawn;
2308 fColorBottom = fColorMain[fEdgeIdx];
2309 }
2310 (this->*fDrawFace)(icodes, xyz, 4, iface, tface);
2311 }
2312 // Draw top face
2313 if (ivis[5] > 0) {
2314 icodes[2] = nv - 1;
2315 icodes[3] = 6;
2316 for (Int_t i = 1; i <= 4; ++i) {
2317 iface[i - 1] = i + 4;
2318 tface[i - 1] = tt[i + (nv << 2) - 5];
2319 }
2320 Int_t cs = fColorTop;
2321 if ( nv <= 3 ) fEdgeIdx = 0 ; // no stack or stack with only one histo
2322 else {
2323 if ( nv > 2 && (v[nv-1] == v[nv-2])) {
2324 for (Int_t iv = nv-1; iv > 2; --iv) {
2325 if (v[nv-1] == v[iv-1]) {
2326 fColorTop = fColorMain[iv-2];
2327 fEdgeIdx = iv - 2;
2328 }
2329 }
2330 }
2331 }
2332 (this->*fDrawFace)(icodes, xyz, 4, iface, tface);
2333 fColorTop = cs;
2334 }
2335 }
2336 }
2337}
2338
2339////////////////////////////////////////////////////////////////////////////////
2340/// Draw stack of lego-plots in polar coordinates
2341///
2342/// \param[in] iordr order of variables (0 - R,PHI; 1 - PHI,R)
2343/// \param[in] na number of steps along 1st variable
2344/// \param[in] nb number of steps along 2nd variable
2345/// \param[in] chopt specific options
2346///
2347/// - `chopt` = 'BF' from BACK to FRONT
2348/// - `chopt` = 'FB' from FRONT to BACK
2349
2350void TPainter3dAlgorithms::LegoPolar(Int_t iordr, Int_t na, Int_t nb, const char *chopt)
2351{
2352
2353 Int_t iphi, jphi, kphi, incr, nphi, ivis[6], iopt, iphi1, iphi2, iface[4], i, j;
2354 Double_t tface[4];
2355 Int_t incrr, k1, k2, ia, ib, ir1, ir2;
2356 Double_t ab[8]; // was [2][4]
2357 Int_t ir, jr, iv, nr, nv, icodes[4];
2358 Double_t xyz[24]; // was [3][8]
2359 ia = ib = 0;
2360 Int_t firstStackNumberDrawn = -1 ; // necessary to compute fColorBottom when the 0 option is set and when the stack is seen from below (bottomview, theta<0.)
2361
2362 TView *view = gPad ? gPad->GetView() : nullptr;
2363 if (!view) {
2364 Error("LegoPolar", "no TView in current pad");
2365 return;
2366 }
2367
2368 if (iordr == 0) {
2369 jr = 1;
2370 jphi = 2;
2371 nr = na;
2372 nphi = nb;
2373 } else {
2374 jr = 2;
2375 jphi = 1;
2376 nr = nb;
2377 nphi = na;
2378 }
2379 if (fNaphi < nphi + 3) {
2380 fNaphi = nphi + 3;
2381 fAphi.resize(fNaphi);
2382 }
2383 if (fAphi.empty()) {
2384 Error("LegoPolar", "failed to allocate array fAphi[%d]", fNaphi);
2385 fNaphi = 0;
2386 return;
2387 }
2388 iopt = 2;
2389 if (*chopt == 'B' || *chopt == 'b') iopt = 1;
2390
2391 // Allocate v and tt arrays
2392 Int_t vSize = fNStack+2;
2393 std::vector<Double_t> v(vSize), tt(4*vSize);
2394
2395 // P R E P A R E P H I A R R A Y
2396 // F I N D C R I T I C A L S E C T O R S
2397 nv = 0;
2398 kphi = nphi;
2399 if (iordr == 0) ia = nr;
2400 if (iordr != 0) ib = nr;
2401 for (i = 1; i <= nphi; ++i) {
2402 if (iordr == 0) ib = i;
2403 if (iordr != 0) ia = i;
2404 (this->*fLegoFunction)(ia, ib, nv, ab, v.data(), tt.data());
2405 if (i == 1) fAphi[0] = ab[jphi - 1];
2406 fAphi[i - 1] = (fAphi[i - 1] + ab[jphi - 1]) / (float)2.;
2407 fAphi[i] = ab[jphi + 3];
2408 }
2409 view->FindPhiSectors(iopt, kphi, fAphi.data(), iphi1, iphi2);
2410
2411 // E N C O D E V I S I B I L I T Y O F S I D E S
2412 // A N D O R D E R A L O N G R
2413 for (i = 1; i <= nphi; ++i) {
2414 if (!iordr) ib = i;
2415 if (iordr) ia = i;
2416 (this->*fLegoFunction)(ia, ib, nv, ab, v.data(), tt.data());
2417 SideVisibilityEncode(iopt, ab[jphi - 1]*kRad, ab[jphi + 3]*kRad, fAphi[i - 1]);
2418 }
2419
2420 // D R A W S T A C K O F L E G O - P L O T S
2421 incr = 1;
2422 iphi = iphi1;
2423L100:
2424 if (iphi > nphi) goto L300;
2425
2426 // D E C O D E V I S I B I L I T Y O F S I D E S
2427 SideVisibilityDecode(fAphi[iphi - 1], ivis[0], ivis[1], ivis[2], ivis[3], ivis[4], ivis[5], incrr);
2428 ir1 = 1;
2429 if (incrr < 0) ir1 = nr;
2430 ir2 = nr - ir1 + 1;
2431 // D R A W L E G O S F O R S E C T O R
2432 for (ir = ir1; incrr < 0 ? ir >= ir2 : ir <= ir2; ir += incrr) {
2433 if (iordr == 0) { ia = ir; ib = iphi; }
2434 else { ia = iphi; ib = ir; }
2435 (this->*fLegoFunction)(ia, ib, nv, ab, v.data(), tt.data());
2436 if (nv < 2 || nv > vSize) continue;
2437 if (Hoption.Zero) {
2438 Double_t total_content=0;
2439 for (iv = 1; iv < nv; ++iv) total_content += v[iv];
2440 if (total_content==0) continue;
2441 }
2442 icodes[0] = ia;
2443 icodes[1] = ib;
2444 for (i = 1; i <= 4; ++i) {
2445 j = i;
2446 if (iordr != 0 && i == 2) j = 4;
2447 if (iordr != 0 && i == 4) j = 2;
2448 xyz[j*3 - 3] = ab[jr + 2*i - 3]*TMath::Cos(ab[jphi + 2*i - 3]*kRad);
2449 xyz[j*3 - 2] = ab[jr + 2*i - 3]*TMath::Sin(ab[jphi + 2*i - 3]*kRad);
2450 xyz[(j + 4)*3 - 3] = xyz[j*3 - 3];
2451 xyz[(j + 4)*3 - 2] = xyz[j*3 - 2];
2452 }
2453 // D R A W S T A C K
2454 firstStackNumberDrawn = -1;
2455 for (iv = 1; iv < nv; ++iv) {
2456 for (i = 1; i <= 4; ++i) {
2457 xyz[i*3 - 1] = v[iv - 1];
2458 xyz[(i + 4)*3 - 1] = v[iv];
2459 }
2460 if (v[iv - 1] >= v[iv]) continue;
2461 icodes[2] = iv;
2462 for (i = 1; i <= 4; ++i) {
2463 if (ivis[i - 1] == 0) continue;
2464 k1 = i - 1;
2465 if (i == 1) k1 = 4;
2466 k2 = i;
2467 if (xyz[k1*3 - 3] == xyz[k2*3 - 3] && xyz[k1*3 - 2] ==
2468 xyz[k2*3 - 2]) continue;
2469 iface[0] = k1;
2470 iface[1] = k2;
2471 iface[2] = k2 + 4;
2472 iface[3] = k1 + 4;
2473 tface[0] = tt[k1 + (iv << 2) - 5];
2474 tface[1] = tt[k2 + (iv << 2) - 5];
2475 tface[2] = tt[k2 + ((iv + 1) << 2) - 5];
2476 tface[3] = tt[k1 + ((iv + 1) << 2) - 5];
2477 icodes[3] = i;
2478 fEdgeIdx = iv-1;
2479 (this->*fDrawFace)(icodes, xyz, 4, iface, tface);
2480 }
2481 if ( firstStackNumberDrawn==-1 ) firstStackNumberDrawn = fEdgeIdx;
2482 }
2483 // D R A W B O T T O M F A C E
2484 if (ivis[4] != 0) {
2485 icodes[2] = 1;
2486 icodes[3] = 5;
2487 for (i = 1; i <= 4; ++i) {
2488 xyz[i*3 - 1] = v[0];
2489 iface[i - 1] = 5 - i;
2490 tface[i - 1] = tt[5 - i - 1];
2491 }
2492 if (!Hoption.Zero) fEdgeIdx = 0;
2493 else {
2494 fEdgeIdx = firstStackNumberDrawn;
2495 fColorBottom = fColorMain[fEdgeIdx];
2496 }
2497 (this->*fDrawFace)(icodes, xyz, 4, iface, tface);
2498 }
2499 // D R A W T O P F A C E
2500 if (ivis[5] != 0) {
2501 icodes[2] = nv - 1;
2502 icodes[3] = 6;
2503 for (i = 1; i <= 4; ++i) {
2504 iface[i - 1] = i + 4;
2505 tface[i - 1] = tt[i + (nv << 2) - 5];
2506 }
2507 Int_t cs = fColorTop;
2508 if ( nv <= 3 ) fEdgeIdx = 0 ; // no stack or stack with only one histo
2509 else {
2510 if ( nv > 2 && (v[nv-1] == v[nv-2])) {
2511 for (iv = nv-1; iv>2; iv--) {
2512 if (v[nv-1] == v[iv-1]) {
2513 fColorTop = fColorMain[iv-2];
2514 fEdgeIdx = iv-2;
2515 }
2516 }
2517 }
2518 }
2519 (this->*fDrawFace)(icodes, xyz, 4, iface, tface);
2520 fColorTop = cs;
2521 }
2522 }
2523 // N E X T P H I
2524L300:
2525 iphi += incr;
2526 if (iphi == 0) iphi = kphi;
2527 if (iphi > kphi) iphi = 1;
2528 if (iphi != iphi2) goto L100;
2529 if (incr == 0)
2530 return;
2531 if (incr < 0) {
2532 incr = 0;
2533 goto L100;
2534 }
2535 incr = -1;
2536 iphi = iphi1;
2537 goto L300;
2538}
2539
2540////////////////////////////////////////////////////////////////////////////////
2541/// Draw stack of lego-plots in cylindrical coordinates
2542///
2543/// \param[in] iordr order of variables (0 - Z,PHI; 1 - PHI,Z)
2544/// \param[in] na number of steps along 1st variable
2545/// \param[in] nb number of steps along 2nd variable
2546/// \param[in] chopt specific options
2547///
2548/// - `chopt` = 'BF' from BACK to FRONT
2549/// - `chopt` = 'FB' from FRONT to BACK
2550
2551void TPainter3dAlgorithms::LegoCylindrical(Int_t iordr, Int_t na, Int_t nb, const char *chopt)
2552{
2553
2554
2555 Int_t iphi, jphi, kphi, incr, nphi, ivis[6], iopt, iphi1, iphi2, iface[4], i, j;
2556 Double_t tface[4], z;
2557 Double_t ab[8]; // was [2][4]
2558 Int_t ia, ib, idummy, iz1, iz2, nz, incrz, k1, k2, nv;
2559 Int_t iv, iz, jz, icodes[4];
2560 Double_t cosphi[4];
2561 Double_t sinphi[4];
2562 Double_t xyz[24]; // was [3][8]
2563 ia = ib = 0;
2564 Int_t firstStackNumberDrawn=-1 ; // necessary to compute fColorBottom when the 0 option is set and when the stack is seen from below (bottomview, theta<0.)
2565
2566 TView *view = gPad ? gPad->GetView() : nullptr;
2567 if (!view) {
2568 Error("LegoCylindrical", "no TView in current pad");
2569 return;
2570 }
2571
2572 if (iordr == 0) {
2573 jz = 1;
2574 jphi = 2;
2575 nz = na;
2576 nphi = nb;
2577 } else {
2578 jz = 2;
2579 jphi = 1;
2580 nz = nb;
2581 nphi = na;
2582 }
2583 if (fNaphi < nphi + 3) {
2584 fNaphi = nphi + 3;
2585 fAphi.resize(fNaphi);
2586 }
2587 if (fAphi.empty()) {
2588 Error("LegoCylindrical", "failed to allocate array fAphi[%d]", fNaphi);
2589 fNaphi = 0;
2590 return;
2591 }
2592 iopt = 2;
2593 if (*chopt == 'B' || *chopt == 'b') iopt = 1;
2594
2595 // Allocate v and tt arrays
2596 Int_t vSize = fNStack+2;
2597 std::vector<Double_t> v(vSize), tt(4*vSize);
2598
2599 // P R E P A R E P H I A R R A Y
2600 // F I N D C R I T I C A L S E C T O R S
2601 nv = 0;
2602 kphi = nphi;
2603 if (iordr == 0) ia = nz;
2604 if (iordr != 0) ib = nz;
2605 for (i = 1; i <= nphi; ++i) {
2606 if (iordr == 0) ib = i;
2607 if (iordr != 0) ia = i;
2608 (this->*fLegoFunction)(ia, ib, nv, ab, v.data(), tt.data());
2609 if (i == 1) fAphi[0] = ab[jphi - 1];
2610 fAphi[i - 1] = (fAphi[i - 1] + ab[jphi - 1]) / (float)2.;
2611 fAphi[i] = ab[jphi + 3];
2612 }
2613 view->FindPhiSectors(iopt, kphi, fAphi.data(), iphi1, iphi2);
2614
2615 // E N C O D E V I S I B I L I T Y O F S I D E S
2616 // A N D O R D E R A L O N G R
2617 for (i = 1; i <= nphi; ++i) {
2618 if (iordr == 0) ib = i;
2619 if (iordr != 0) ia = i;
2620 (this->*fLegoFunction)(ia, ib, nv, ab, v.data(), tt.data());
2621 SideVisibilityEncode(iopt, ab[jphi - 1]*kRad, ab[jphi + 3]*kRad, fAphi[i - 1]);
2622 }
2623
2624 // F I N D O R D E R A L O N G Z
2625 incrz = 1;
2626 iz1 = 1;
2627 view->FindNormal(0, 0, 1, z);
2628 if ((z <= 0 && iopt == 1) || (z > 0 && iopt == 2)) {
2629 incrz = -1;
2630 iz1 = nz;
2631 }
2632 iz2 = nz - iz1 + 1;
2633
2634 // D R A W S T A C K O F L E G O - P L O T S
2635 incr = 1;
2636 iphi = iphi1;
2637L100:
2638 if (iphi > nphi) goto L400;
2639 // D E C O D E V I S I B I L I T Y O F S I D E S
2640 idummy = 0;
2641 SideVisibilityDecode(fAphi[iphi - 1], ivis[4], ivis[1], ivis[5], ivis[3], ivis[0], ivis[2], idummy);
2642 for (iz = iz1; incrz < 0 ? iz >= iz2 : iz <= iz2; iz += incrz) {
2643 if (iordr == 0) {ia = iz; ib = iphi;}
2644 else {ia = iphi; ib = iz;}
2645 (this->*fLegoFunction)(ia, ib, nv, ab, v.data(), tt.data());
2646 if (nv < 2 || nv > vSize) continue;
2647 icodes[0] = ia;
2648 icodes[1] = ib;
2649 for (i = 1; i <= 4; ++i) {
2650 j = i;
2651 if (iordr != 0 && i == 2) j = 4;
2652 if (iordr != 0 && i == 4) j = 2;
2653 cosphi[j - 1] = TMath::Cos(ab[jphi + 2*i - 3]*kRad);
2654 sinphi[j - 1] = TMath::Sin(ab[jphi + 2*i - 3]*kRad);
2655 xyz[j*3 - 1] = ab[jz + 2*i - 3];
2656 xyz[(j + 4)*3 - 1] = ab[jz + 2*i - 3];
2657 }
2658 // D R A W S T A C K
2659 firstStackNumberDrawn = -1;
2660 for (iv = 1; iv < nv; ++iv) {
2661 for (i = 1; i <= 4; ++i) {
2662 xyz[i*3 - 3] = v[iv - 1]*cosphi[i - 1];
2663 xyz[i*3 - 2] = v[iv - 1]*sinphi[i - 1];
2664 xyz[(i + 4)*3 - 3] = v[iv]*cosphi[i - 1];
2665 xyz[(i + 4)*3 - 2] = v[iv]*sinphi[i - 1];
2666 }
2667 if (v[iv - 1] >= v[iv]) continue;
2668 icodes[2] = iv;
2669 for (i = 1; i <= 4; ++i) {
2670 if (ivis[i - 1] == 0) continue;
2671 k1 = i;
2672 k2 = i - 1;
2673 if (i == 1) k2 = 4;
2674 iface[0] = k1;
2675 iface[1] = k2;
2676 iface[2] = k2 + 4;
2677 iface[3] = k1 + 4;
2678 tface[0] = tt[k1 + (iv << 2) - 5];
2679 tface[1] = tt[k2 + (iv << 2) - 5];
2680 tface[2] = tt[k2 + ((iv + 1) << 2) - 5];
2681 tface[3] = tt[k1 + ((iv + 1) << 2) - 5];
2682 icodes[3] = i;
2683 fEdgeIdx = iv-1;
2684 (this->*fDrawFace)(icodes, xyz, 4, iface, tface);
2685 }
2686 if ( firstStackNumberDrawn==-1 ) firstStackNumberDrawn = fEdgeIdx;
2687 }
2688 // D R A W B O T T O M F A C E
2689 if (ivis[4] != 0 && v[0] > 0) {
2690 icodes[2] = 1;
2691 icodes[3] = 5;
2692 for (i = 1; i <= 4; ++i) {
2693 xyz[i*3 - 3] = v[0]*cosphi[i - 1];
2694 xyz[i*3 - 2] = v[0]*sinphi[i - 1];
2695 iface[i - 1] = i;
2696 tface[i - 1] = tt[i - 1];
2697 }
2698 if (!Hoption.Zero) fEdgeIdx = 0;
2699 else {
2700 fEdgeIdx = firstStackNumberDrawn;
2701 fColorBottom = fColorMain[fEdgeIdx];
2702 }
2703 (this->*fDrawFace)(icodes, xyz, 4, iface, tface);
2704 }
2705 // D R A W T O P F A C E
2706 if (ivis[5] != 0 && v[nv - 1] > 0) {
2707 icodes[2] = nv - 1;
2708 icodes[3] = 6;
2709 for (i = 1; i <= 4; ++i) {
2710 iface[i - 1] = 5 - i + 4;
2711 tface[i - 1] = tt[5 - i + (nv << 2) - 5];
2712 }
2713 Int_t cs = fColorTop;
2714 if ( nv <= 3 ) fEdgeIdx = 0 ; // no stack or stack with only one histo
2715 else {
2716 if ( nv > 2 && (v[nv-1] == v[nv-2])) {
2717 for (iv = nv-1; iv>2; iv--) {
2718 if (v[nv-1] == v[iv-1]) {
2719 fColorTop = fColorMain[iv-2];
2720 fEdgeIdx = iv-2;
2721 }
2722 }
2723 }
2724 }
2725 (this->*fDrawFace)(icodes, xyz, 4, iface, tface);
2726 fColorTop = cs;
2727 }
2728 }
2729 // N E X T P H I
2730L400:
2731 iphi += incr;
2732 if (iphi == 0) iphi = kphi;
2733 if (iphi > kphi) iphi = 1;
2734 if (iphi != iphi2) goto L100;
2735 if (incr == 0)
2736 return;
2737 if (incr < 0) {
2738 incr = 0;
2739 goto L100;
2740 }
2741 incr = -1;
2742 iphi = iphi1;
2743 goto L400;
2744}
2745
2746////////////////////////////////////////////////////////////////////////////////
2747/// Draw stack of lego-plots spheric coordinates
2748///
2749/// \param[in] ipsdr pseudo-rapidity flag
2750/// \param[in] iordr order of variables (0 - THETA,PHI; 1 - PHI,THETA)
2751/// \param[in] na number of steps along 1st variable
2752/// \param[in] nb number of steps along 2nd variable
2753/// \param[in] chopt specific options
2754///
2755/// - `chopt` = 'BF' from BACK to FRONT
2756/// - `chopt` = 'FB' from FRONT to BACK
2757
2758void TPainter3dAlgorithms::LegoSpherical(Int_t ipsdr, Int_t iordr, Int_t na, Int_t nb, const char *chopt)
2759{
2760 Int_t iphi, jphi, kphi, incr, nphi, ivis[6], iopt, iphi1, iphi2, iface[4], i, j;
2761 Double_t tface[4], costh[4];
2762 Double_t sinth[4];
2763 Int_t k1, k2, ia, ib, incrth, ith, jth, kth, nth, mth, ith1, ith2, nv;
2764 Double_t ab[8]; // was [2][4]
2765 Double_t th;
2766 Int_t iv, icodes[4];
2767 Double_t zn, cosphi[4];
2768 Double_t sinphi[4], th1, th2, phi;
2769 Double_t xyz[24]; // was [3][8]
2770 Double_t phi1, phi2;
2771 ia = ib = 0;
2772 Int_t firstStackNumberDrawn=-1 ; // necessary to compute fColorBottom when the 0 option is set and when the stack is seen from below (bottomview, theta<0.)
2773
2774 TView *view = gPad ? gPad->GetView() : nullptr;
2775 if (!view) {
2776 Error("LegoSpherical", "no TView in current pad");
2777 return;
2778 }
2779
2780 if (iordr == 0) {
2781 jth = 1;
2782 jphi = 2;
2783 nth = na;
2784 nphi = nb;
2785 } else {
2786 jth = 2;
2787 jphi = 1;
2788 nth = nb;
2789 nphi = na;
2790 }
2791 if (fNaphi < nth + 3 || fNaphi < nphi + 3) {
2792 fNaphi = TMath::Max(nth, nphi) + 3;
2793 fAphi.resize(fNaphi);
2794 }
2795 if (fAphi.empty()) {
2796 Error("LegoSpherical", "failed to allocate array fAphi[%d]", fNaphi);
2797 fNaphi = 0;
2798 return;
2799 }
2800 iopt = 2;
2801 if (*chopt == 'B' || *chopt == 'b') iopt = 1;
2802
2803 // Allocate v and tt arrays
2804 Int_t vSize = fNStack+2;
2805 std::vector<Double_t> v(vSize), tt(4*vSize);
2806
2807 // P R E P A R E P H I A R R A Y
2808 // F I N D C R I T I C A L P H I S E C T O R S
2809 nv = 0;
2810 kphi = nphi;
2811 mth = nth / 2;
2812 if (mth == 0) mth = 1;
2813 if (iordr == 0) ia = mth;
2814 if (iordr != 0) ib = mth;
2815 for (i = 1; i <= nphi; ++i) {
2816 if (iordr == 0) ib = i;
2817 if (iordr != 0) ia = i;
2818 (this->*fLegoFunction)(ia, ib, nv, ab, v.data(), tt.data());
2819 if (i == 1) fAphi[0] = ab[jphi - 1];
2820 fAphi[i - 1] = (fAphi[i - 1] + ab[jphi - 1]) / (float)2.;
2821 fAphi[i] = ab[jphi + 3];
2822 }
2823 view->FindPhiSectors(iopt, kphi, fAphi.data(), iphi1, iphi2);
2824
2825 // P R E P A R E T H E T A A R R A Y
2826 if (iordr == 0) ib = 1;
2827 if (iordr != 0) ia = 1;
2828 for (i = 1; i <= nth; ++i) {
2829 if (iordr == 0) ia = i;
2830 if (iordr != 0) ib = i;
2831 (this->*fLegoFunction)(ia, ib, nv, ab, v.data(), tt.data());
2832 if (i == 1) fAphi[0] = ab[jth - 1];
2833 fAphi[i - 1] = (fAphi[i - 1] + ab[jth - 1]) / (float)2.;
2834 fAphi[i] = ab[jth + 3];
2835 }
2836
2837 // D R A W S T A C K O F L E G O - P L O T S
2838 kth = nth;
2839
2840 incr = 1;
2841 iphi = iphi1;
2842L100:
2843 if (iphi > nphi) goto L500;
2844
2845 // F I N D C R I T I C A L T H E T A S E C T O R S
2846 if (!iordr) {ia = mth; ib = iphi; }
2847 else {ia = iphi; ib = mth; }
2848 (this->*fLegoFunction)(ia, ib, nv, ab, v.data(), tt.data());
2849 phi = (ab[jphi - 1] + ab[jphi + 3]) / (float)2.;
2850 view->FindThetaSectors(iopt, phi, kth, fAphi.data(), ith1, ith2);
2851 incrth = 1;
2852 ith = ith1;
2853L200:
2854 if (ith > nth) goto L400;
2855 if (iordr == 0) ia = ith;
2856 if (iordr != 0) ib = ith;
2857 (this->*fLegoFunction)(ia, ib, nv, ab, v.data(), tt.data());
2858 if (nv < 2 || nv > vSize) goto L400;
2859
2860 // D E F I N E V I S I B I L I T Y O F S I D E S
2861 for (i = 1; i <= 6; ++i) ivis[i - 1] = 0;
2862
2863 phi1 = kRad*ab[jphi - 1];
2864 phi2 = kRad*ab[jphi + 3];
2865 th1 = kRad*ab[jth - 1];
2866 th2 = kRad*ab[jth + 3];
2867 view->FindNormal(TMath::Sin(phi1), -TMath::Cos(phi1), 0, zn);
2868 if (zn > 0) ivis[1] = 1;
2869 view->FindNormal(-TMath::Sin(phi2), TMath::Cos(phi2), 0, zn);
2870 if (zn > 0) ivis[3] = 1;
2871 phi = (phi1 + phi2) / (float)2.;
2872 view->FindNormal(-TMath::Cos(phi)*TMath::Cos(th1), -TMath::Sin(phi)*TMath::Cos(th1), TMath::Sin(th1), zn);
2873 if (zn > 0) ivis[0] = 1;
2874 view->FindNormal(TMath::Cos(phi)*TMath::Cos(th2), TMath::Sin(phi)*TMath::Cos(th2), -TMath::Sin(th2), zn);
2875 if (zn > 0) ivis[2] = 1;
2876 th = (th1 + th2) / (float)2.;
2877 if (ipsdr == 1) th = kRad*90;
2878 view->FindNormal(TMath::Cos(phi)*TMath::Sin(th), TMath::Sin(phi)*TMath::Sin(th), TMath::Cos(th), zn);
2879 if (zn < 0) ivis[4] = 1;
2880 if (zn > 0) ivis[5] = 1;
2881
2882 // D R A W S T A C K
2883 icodes[0] = ia;
2884 icodes[1] = ib;
2885 for (i = 1; i <= 4; ++i) {
2886 j = i;
2887 if (iordr != 0 && i == 2) j = 4;
2888 if (iordr != 0 && i == 4) j = 2;
2889 costh[j - 1] = TMath::Cos(kRad*ab[jth + 2*i - 3]);
2890 sinth[j - 1] = TMath::Sin(kRad*ab[jth + 2*i - 3]);
2891 cosphi[j - 1] = TMath::Cos(kRad*ab[jphi + 2*i - 3]);
2892 sinphi[j - 1] = TMath::Sin(kRad*ab[jphi + 2*i - 3]);
2893 }
2894 firstStackNumberDrawn = -1;
2895 for (iv = 1; iv < nv; ++iv) {
2896 if (ipsdr == 1) {
2897 for (i = 1; i <= 4; ++i) {
2898 xyz[i*3 - 3] = v[iv - 1]*cosphi[i - 1];
2899 xyz[i*3 - 2] = v[iv - 1]*sinphi[i - 1];
2900 xyz[i*3 - 1] = v[iv - 1]*costh[i - 1] / sinth[i - 1];
2901 xyz[(i + 4)*3 - 3] = v[iv]*cosphi[i - 1];
2902 xyz[(i + 4)*3 - 2] = v[iv]*sinphi[i - 1];
2903 xyz[(i + 4)*3 - 1] = v[iv]*costh[i - 1] / sinth[i - 1];
2904 }
2905 } else {
2906 for (i = 1; i <= 4; ++i) {
2907 xyz[i*3 - 3] = v[iv - 1]*sinth[i - 1]*cosphi[i - 1];
2908 xyz[i*3 - 2] = v[iv - 1]*sinth[i - 1]*sinphi[i - 1];
2909 xyz[i*3 - 1] = v[iv - 1]*costh[i - 1];
2910 xyz[(i + 4)*3 - 3] = v[iv]*sinth[i - 1]*cosphi[i - 1];
2911 xyz[(i + 4)*3 - 2] = v[iv]*sinth[i - 1]*sinphi[i - 1];
2912 xyz[(i + 4)*3 - 1] = v[iv]*costh[i - 1];
2913 }
2914 }
2915 if (v[iv - 1] >= v[iv]) continue;
2916 icodes[2] = iv;
2917 for (i = 1; i <= 4; ++i) {
2918 if (ivis[i - 1] == 0) continue;
2919 k1 = i - 1;
2920 if (i == 1) k1 = 4;
2921 k2 = i;
2922 iface[0] = k1;
2923 iface[1] = k2;
2924 iface[2] = k2 + 4;
2925 iface[3] = k1 + 4;
2926 tface[0] = tt[k1 + (iv << 2) - 5];
2927 tface[1] = tt[k2 + (iv << 2) - 5];
2928 tface[2] = tt[k2 + ((iv + 1) << 2) - 5];
2929 tface[3] = tt[k1 + ((iv + 1) << 2) - 5];
2930 icodes[3] = i;
2931 fEdgeIdx = iv-1;
2932 (this->*fDrawFace)(icodes, xyz, 4, iface, tface);
2933 }
2934 if ( firstStackNumberDrawn==-1 ) firstStackNumberDrawn = fEdgeIdx;
2935 }
2936 // D R A W B O T T O M F A C E
2937 if (ivis[4] != 0 && v[0] > 0) {
2938 icodes[2] = 1;
2939 icodes[3] = 5;
2940 for (i = 1; i <= 4; ++i) {
2941 if (ipsdr == 1) {
2942 xyz[i*3 - 3] = v[0]*cosphi[i - 1];
2943 xyz[i*3 - 2] = v[0]*sinphi[i - 1];
2944 xyz[i*3 - 1] = v[0]*costh[i - 1] / sinth[i - 1];
2945 } else {
2946 xyz[i*3 - 3] = v[0]*sinth[i - 1]*cosphi[i - 1];
2947 xyz[i*3 - 2] = v[0]*sinth[i - 1]*sinphi[i - 1];
2948 xyz[i*3 - 1] = v[0]*costh[i - 1];
2949 }
2950 iface[i - 1] = 5 - i;
2951 tface[i - 1] = tt[5 - i - 1];
2952 }
2953 if (!Hoption.Zero) fEdgeIdx = 0;
2954 else {
2955 fEdgeIdx = firstStackNumberDrawn;
2956 fColorBottom = fColorMain[fEdgeIdx];
2957 }
2958 (this->*fDrawFace)(icodes, xyz, 4, iface, tface);
2959 }
2960 // D R A W T O P F A C E
2961 if (ivis[5] != 0 && v[nv - 1] > 0) {
2962 icodes[2] = nv - 1;
2963 icodes[3] = 6;
2964 for (i = 1; i <= 4; ++i) {
2965 iface[i - 1] = i + 4;
2966 tface[i - 1] = tt[i + 4 + 2*nv - 5];
2967 }
2968 Int_t cs = fColorTop;
2969 if ( nv <= 3 ) fEdgeIdx = 0 ; // no stack or stack with only one histo
2970 else {
2971 if ( nv > 2 && (v[nv-1] == v[nv-2])) {
2972 for (iv = nv-1; iv>2; iv--) {
2973 if (v[nv-1] == v[iv-1]) {
2974 fColorTop = fColorMain[iv-2];
2975 fEdgeIdx = iv-2;
2976 }
2977 }
2978 }
2979 }
2980 (this->*fDrawFace)(icodes, xyz, 4, iface, tface);
2981 fColorTop = cs;
2982 }
2983 // N E X T T H E T A
2984L400:
2985 ith += incrth;
2986 if (ith == 0) ith = kth;
2987 if (ith > kth) ith = 1;
2988 if (ith != ith2) goto L200;
2989 if (incrth == 0) goto L500;
2990 if (incrth < 0) {
2991 incrth = 0;
2992 goto L200;
2993 }
2994 incrth = -1;
2995 ith = ith1;
2996 goto L400;
2997 // N E X T P H I
2998L500:
2999 iphi += incr;
3000 if (iphi == 0) iphi = kphi;
3001 if (iphi > kphi) iphi = 1;
3002 if (iphi != iphi2) goto L100;
3003 if (incr == 0)
3004 return;
3005 if (incr < 0) {
3006 incr = 0;
3007 goto L100;
3008 }
3009 incr = -1;
3010 iphi = iphi1;
3011 goto L500;
3012}
3013
3014////////////////////////////////////////////////////////////////////////////////
3015/// Set light source
3016///
3017/// \param[in] nl source number: 1 off all light sources, 0 set diffused light
3018/// \param[in] yl intensity of the light source
3019/// \param[in] xscr, yscr, zscr direction of the light (in respect of the screen)
3020///
3021/// \param[out] irep reply (0 - O.K, -1 error)
3022
3023
3024void TPainter3dAlgorithms::LightSource(Int_t nl, Double_t yl, Double_t xscr,
3025 Double_t yscr, Double_t zscr, Int_t &irep)
3026{
3027 /* Local variables */
3028 Int_t i;
3029 Double_t s;
3030
3031 irep = 0;
3032 if (nl < 0) goto L100;
3033 else if (nl == 0) goto L200;
3034 else goto L300;
3035
3036 // S W I T C H O F F L I G H T S
3037L100:
3038 fLoff = 1;
3039 fYdl = 0;
3040 for (i = 1; i <= 4; ++i) {
3041 fYls[i - 1] = 0;
3042 }
3043 return;
3044 // S E T D I F F U S E D L I G H T
3045L200:
3046 if (yl < 0) {
3047 Error("LightSource", "negative light intensity");
3048 irep = -1;
3049 return;
3050 }
3051 fYdl = yl;
3052 goto L400;
3053 // S E T L I G H T S O U R C E
3054L300:
3055 if (nl > 4 || yl < 0) {
3056 Error("LightSource", "illegal light source number (nl=%d, yl=%f)", nl, yl);
3057 irep = -1;
3058 return;
3059 }
3060 s = TMath::Sqrt(xscr*xscr + yscr*yscr + zscr*zscr);
3061 if (s == 0) {
3062 Error("LightSource", "light source is placed at origin");
3063 irep = -1;
3064 return;
3065 }
3066 fYls[nl - 1] = yl;
3067 fVls[nl*3 - 3] = xscr / s;
3068 fVls[nl*3 - 2] = yscr / s;
3069 fVls[nl*3 - 1] = zscr / s;
3070 // C H E C K L I G H T S
3071L400:
3072 fLoff = 0;
3073 if (fYdl != 0) return;
3074 for (i = 1; i <= 4; ++i) {
3075 if (fYls[i - 1] != 0) return;
3076 }
3077 fLoff = 1;
3078}
3079
3080////////////////////////////////////////////////////////////////////////////////
3081/// Find surface luminosity at given point
3082///
3083/// \param[in] view pointer on TView object
3084/// \param[in] anorm surface normal at given point
3085///
3086/// \param[out] flum luminosity
3087
3088void TPainter3dAlgorithms::Luminosity(TView *view, Double_t *anorm, Double_t &flum)
3089{
3090 flum = 0;
3091
3092 if (!view || fLoff) return;
3093
3094 /* Local variables */
3095 Double_t cosn, cosr;
3096 Int_t i;
3097 Double_t s, vl[3], vn[3];
3098
3099 // T R A N S F E R N O R M A L T O SCREEN COORDINATES
3100 view->NormalWCtoNDC(anorm, vn);
3101 s = TMath::Sqrt(vn[0]*vn[0] + vn[1]*vn[1] + vn[2]*vn[2]);
3102 if (vn[2] < 0) s = -(Double_t)s;
3103 vn[0] /= s;
3104 vn[1] /= s;
3105 vn[2] /= s;
3106
3107 // F I N D L U M I N O S I T Y
3108 flum = fYdl*fQA;
3109 for (i = 1; i <= 4; ++i) {
3110 if (fYls[i - 1] <= 0) continue;
3111 vl[0] = fVls[i*3 - 3];
3112 vl[1] = fVls[i*3 - 2];
3113 vl[2] = fVls[i*3 - 1];
3114 cosn = vl[0]*vn[0] + vl[1]*vn[1] + vl[2]*vn[2];
3115 if (cosn < 0) continue;
3116 cosr = vn[1]*(vn[2]*vl[1] - vn[1]*vl[2]) - vn[0]*(vn[0]*vl[2]
3117 - vn[2]*vl[0]) + vn[2]*cosn;
3118 if (cosr <= 0) cosr = 0;
3119 flum += fYls[i - 1]*(fQD*cosn + fQS*TMath::Power(cosr, fNqs));
3120 }
3121}
3122
3123////////////////////////////////////////////////////////////////////////////////
3124/// Modify SCREEN
3125///
3126/// \param[in] r1 1-st point of the line
3127/// \param[in] r2 2-nd point of the line
3128
3129void TPainter3dAlgorithms::ModifyScreen(Double_t *r1, Double_t *r2)
3130{
3131 /* Local variables */
3132 Int_t i, i1, i2;
3133 Double_t x1, x2, y1, y2, dy, ww, yy1, yy2, *tn;
3134
3135 /* Parameter adjustments */
3136 --r2;
3137 --r1;
3138
3139 TView *view = gPad ? gPad->GetView() : nullptr;
3140
3141 if (view) {
3142 tn = view->GetTN();
3143 if (tn) {
3144 x1 = tn[0]*r1[1] + tn[1]*r1[2] + tn[2]*r1[3] + tn[3];
3145 x2 = tn[0]*r2[1] + tn[1]*r2[2] + tn[2]*r2[3] + tn[3];
3146 y1 = tn[4]*r1[1] + tn[5]*r1[2] + tn[6]*r1[3] + tn[7];
3147 y2 = tn[4]*r2[1] + tn[5]*r2[2] + tn[6]*r2[3] + tn[7];
3148 } else {
3149 Error("ModifyScreen", "invalid TView in current pad");
3150 return;
3151 }
3152 } else {
3153 Error("ModifyScreen", "no TView in current pad");
3154 return;
3155 }
3156
3157 if (x1 >= x2) {
3158 ww = x1;
3159 x1 = x2;
3160 x2 = ww;
3161 ww = y1;
3162 y1 = y2;
3163 y2 = ww;
3164 }
3165 i1 = Int_t((x1 - fX0) / fDX) + 15;
3166 i2 = Int_t((x2 - fX0) / fDX) + 15;
3167 if (i1 == i2) return;
3168
3169 // M O D I F Y B O U N D A R I E S OF THE SCREEN
3170 dy = (y2 - y1) / (i2 - i1);
3171 for (i = i1; i <= i2 - 1; ++i) {
3172 yy1 = y1 + dy*(i - i1);
3173 yy2 = yy1 + dy;
3174 if (fD[2*i - 2] > yy1) fD[2*i - 2] = yy1;
3175 if (fD[2*i - 1] > yy2) fD[2*i - 1] = yy2;
3176 if (fU[2*i - 2] < yy1) fU[2*i - 2] = yy1;
3177 if (fU[2*i - 1] < yy2) fU[2*i - 1] = yy2;
3178 }
3179}
3180
3181////////////////////////////////////////////////////////////////////////////////
3182/// Store pointer to current algorithm to draw faces
3183
3184void TPainter3dAlgorithms::SetDrawFace(DrawFaceFunc_t drface)
3185{
3186 fDrawFace = drface;
3187}
3188
3189////////////////////////////////////////////////////////////////////////////////
3190/// Store pointer to current lego function
3191
3192void TPainter3dAlgorithms::SetLegoFunction(LegoFunc_t fun)
3193{
3194 fLegoFunction = fun;
3195}
3196
3197////////////////////////////////////////////////////////////////////////////////
3198/// Store pointer to current surface function
3199
3204
3205////////////////////////////////////////////////////////////////////////////////
3206/// Store dark color for stack number n
3207
3208void TPainter3dAlgorithms::SetColorDark(Color_t color, Int_t n)
3209{
3210 if (n < 0 ) {fColorBottom = color; return;}
3211 if (n > fNStack ) {fColorTop = color; return;}
3212 fColorDark[n] = color;
3213}
3214
3215////////////////////////////////////////////////////////////////////////////////
3216/// Store color for stack number n
3217
3218void TPainter3dAlgorithms::SetColorMain(Color_t color, Int_t n)
3219{
3220 if (n < 0 ) {fColorBottom = color; return;}
3221 if (n > fNStack ) {fColorTop = color; return;}
3222 fColorMain[n] = color;
3223}
3224
3225////////////////////////////////////////////////////////////////////////////////
3226
3227void TPainter3dAlgorithms::SetEdgeAtt(Color_t color, Style_t style, Width_t width, Int_t n)
3228{
3229 // Store edge attributes
3230
3231 fEdgeColor[n] = color;
3232 fEdgeStyle[n] = style;
3233 fEdgeWidth[n] = width;
3234}
3235
3236////////////////////////////////////////////////////////////////////////////////
3237/// Decode side visibilities and order along R for sector
3238///
3239/// \param[in] val encoded value
3240///
3241/// \param[out] iv1,iv2,iv3,iv4,iv5,iv6 visibility of the sides
3242/// \param[out] ir increment along R
3243
3244void TPainter3dAlgorithms::SideVisibilityDecode(Double_t val, Int_t &iv1, Int_t &iv2, Int_t &iv3, Int_t &iv4, Int_t &iv5, Int_t &iv6, Int_t &ir)
3245{
3246 Int_t ivis[6], i, k, num;
3247
3248 k = Int_t(val);
3249 num = 128;
3250 for (i = 1; i <= 6; ++i) {
3251 ivis[i - 1] = 0;
3252 num /= 2;
3253 if (k < num) continue;
3254 k -= num;
3255 ivis[i - 1] = 1;
3256 }
3257 ir = 1;
3258 if (k == 1) ir = -1;
3259 iv1 = ivis[5];
3260 iv2 = ivis[4];
3261 iv3 = ivis[3];
3262 iv4 = ivis[2];
3263 iv5 = ivis[1];
3264 iv6 = ivis[0];
3265}
3266
3267////////////////////////////////////////////////////////////////////////////////
3268/// Encode side visibilities and order along R for sector
3269///
3270/// \param[in] iopt options: 1: from BACK to FRONT 'BF', 2: from FRONT to BACK 'FB'
3271/// \param[in] phi1 1st phi of sector
3272/// \param[in] phi2 2nd phi of sector
3273///
3274/// \param[out] val encoded value
3275
3276void TPainter3dAlgorithms::SideVisibilityEncode(Int_t iopt, Double_t phi1, Double_t phi2, Double_t &val)
3277{
3278 /* Local variables */
3279 Double_t zn, phi;
3280 Int_t k = 0;
3281
3282 TView *view = gPad ? gPad->GetView() : nullptr;
3283 if (!view) {
3284 Error("SideVisibilityEncode", "no TView in current pad");
3285 return;
3286 }
3287
3288 view->FindNormal(0, 0, 1, zn);
3289 if (zn > 0) k += 64;
3290 if (zn < 0) k += 32;
3291 view->FindNormal(-TMath::Sin(phi2), TMath::Cos(phi2), 0, zn);
3292 if (zn > 0) k += 16;
3293 view->FindNormal(TMath::Sin(phi1), -TMath::Cos(phi1), 0, zn);
3294 if (zn > 0) k += 4;
3295 phi = (phi1 + phi2) / (float)2.;
3296 view->FindNormal(TMath::Cos(phi), TMath::Sin(phi), 0, zn);
3297 if (zn > 0) k += 8;
3298 if (zn < 0) k += 2;
3299 if ((zn <= 0 && iopt == 1) || (zn > 0 && iopt == 2)) ++k;
3300 val = Double_t(k);
3301}
3302
3303////////////////////////////////////////////////////////////////////////////////
3304/// Set Spectrum
3305///
3306/// \param[in] nl number of levels
3307/// \param[in] fmin MIN function value
3308/// \param[in] fmax MAX function value
3309/// \param[in] ic initial color index (for 1st level)
3310/// \param[in] idc color index increment
3311///
3312/// \param[out] irep reply (0 O.K., -1 error)
3313
3314void TPainter3dAlgorithms::Spectrum(Int_t nl, Double_t fmin, Double_t fmax, Int_t ic, Int_t idc, Int_t &irep)
3315{
3316 static const char *where = "Spectrum";
3317
3318 /* Local variables */
3319 Double_t delf;
3320 Int_t i;
3321
3322 irep = 0;
3323 if (nl == 0) {fNlevel = 0; return; }
3324
3325 // C H E C K P A R A M E T E R S
3326 if (fmax <= fmin) {
3327 Error(where, "fmax (%f) less than fmin (%f)", fmax, fmin);
3328 irep = -1;
3329 return;
3330 }
3331 if (nl < 0 || nl > 256) {
3332 Error(where, "illegal number of levels (%d)", nl);
3333 irep = -1;
3334 return;
3335 }
3336 if (ic < 0) {
3337 Error(where, "initial color index is negative");
3338 irep = -1;
3339 return;
3340 }
3341 if (idc < 0) {
3342 Error(where, "color index increment must be positive");
3343 irep = -1;
3344 }
3345
3346 // S E T S P E C T R
3347 const Int_t kMAXCOL = 50;
3348 delf = (fmax - fmin) / nl;
3349 fNlevel = -(nl + 1);
3350 for (i = 1; i <= nl+1; ++i) {
3351 fFunLevel[i - 1] = fmin + (i - 1)*delf;
3352 fColorLevel[i] = ic + (i - 1)*idc;
3353 if (ic <= kMAXCOL && fColorLevel[i] > kMAXCOL) fColorLevel[i] -= kMAXCOL;
3354 }
3355 fColorLevel[0] = fColorLevel[1];
3356 fColorLevel[nl + 1] = fColorLevel[nl];
3357}
3358
3359////////////////////////////////////////////////////////////////////////////////
3360/// Draw surface in cartesian coordinate system
3361///
3362/// \param[in] ang angle between X ang Y (not used in this method)
3363/// \param[in] nx number of steps along X
3364/// \param[in] ny number of steps along Y
3365/// \param[in] chopt specific options
3366///
3367/// - `chopt` = 'BF' from BACK to FRONT
3368/// - `chopt` = 'FB' from FRONT to BACK
3369
3370void TPainter3dAlgorithms::SurfaceCartesian(Double_t, Int_t nx, Int_t ny, const char *chopt)
3371{
3372 Int_t iface[4] = { 1,2,3,4 };
3373 Int_t icodes[3];
3374 Double_t f[4*3], tt[4], xyz[4*3];
3375
3376 TView *view = gPad ? gPad->GetView() : nullptr;
3377 if (!view) {
3378 Error("SurfaceCartesian", "no TView in current pad");
3379 return;
3380 }
3381 Double_t *tnorm = view->GetTnorm();
3382 if (!tnorm) return;
3383
3384 // Define order of drawing
3385 Int_t incrx = (tnorm[8] < 0.) ? -1 : +1;
3386 Int_t incry = (tnorm[9] < 0.) ? -1 : +1;
3387 if (*chopt != 'B' && *chopt != 'b') { // front to back
3388 incrx = -incrx; incry = -incry;
3389 }
3390 Int_t ix1 = (incrx == +1) ? 1 : nx;
3391 Int_t iy1 = (incry == +1) ? 1 : ny;
3392 Int_t ix2 = (incrx == +1) ? nx : 1;
3393 Int_t iy2 = (incry == +1) ? ny : 1;
3394
3395 // Draw surface
3396 THistPainter *painter = (THistPainter*)gCurrentHist->GetPainter();
3397 for (Int_t iy = iy1; iy != iy2+incry; iy += incry) {
3398 for (Int_t ix = ix1; ix != ix2+incrx; ix += incrx) {
3399 if (!painter->IsInside(ix,iy)) continue;
3400 (this->*fSurfaceFunction)(ix, iy, f, tt);
3401 for (Int_t i = 0; i < 4; ++i) {
3402 xyz[i*3 + 0] = f[i*3 + 0];
3403 xyz[i*3 + 1] = f[i*3 + 1];
3404 xyz[i*3 + 2] = f[i*3 + 2];
3405 // added EJB -->
3406 Double_t al, ab;
3407 if (Hoption.Proj == 1 ) {
3408 THistPainter::ProjectAitoff2xy(xyz[i*3 + 0], xyz[i*3 + 1], al, ab);
3409 xyz[i*3 + 0] = al;
3410 xyz[i*3 + 1] = ab;
3411 } else if (Hoption.Proj == 2 ) {
3412 THistPainter::ProjectMercator2xy(xyz[i*3 + 0], xyz[i*3 + 1], al, ab);
3413 xyz[i*3 + 0] = al;
3414 xyz[i*3 + 1] = ab;
3415 } else if (Hoption.Proj == 3) {
3416 THistPainter::ProjectSinusoidal2xy(xyz[i*3 + 0], xyz[i*3 + 1], al, ab);
3417 xyz[i*3 + 0] = al;
3418 xyz[i*3 + 1] = ab;
3419 } else if (Hoption.Proj == 4) {
3420 THistPainter::ProjectParabolic2xy(xyz[i*3 + 0], xyz[i*3 + 1], al, ab);
3421 xyz[i*3 + 0] = al;
3422 xyz[i*3 + 1] = ab;
3423 } else if (Hoption.Proj == 5) {
3424 THistPainter::ProjectMollweide2xy(xyz[i*3 + 0], xyz[i*3 + 1], al, ab);
3425 xyz[i*3 + 0] = al;
3426 xyz[i*3 + 1] = ab;
3427 }
3428 }
3429 icodes[0] = ix;
3430 icodes[1] = iy;
3431 icodes[2] = -1; // -1 for data, 0 for front a back boxes
3432 fEdgeIdx = 0; // constant since stacks are not (yet?) handled for surfaces
3433 (this->*fDrawFace)(icodes, xyz, 4, iface, tt);
3434 }
3435 }
3436}
3437
3438////////////////////////////////////////////////////////////////////////////////
3439/// Service function for Surfaces
3440
3441void TPainter3dAlgorithms::SurfaceFunction(Int_t ia, Int_t ib, Double_t *f, Double_t *t)
3442{
3443 static Int_t ixadd[4] = { 0,1,1,0 };
3444 static Int_t iyadd[4] = { 0,0,1,1 };
3445
3446 Double_t rinrad = gStyle->GetLegoInnerR();
3447 Double_t dangle = 10; //Delta angle for Rapidity option
3448 Double_t yval1l, yval2l;
3449 Double_t xlab1l, xlab2l, ylab1l, ylab2l;
3450 Int_t i, ixa, iya, icx, ixt, iyt;
3451
3452 /* Parameter adjustments */
3453 --t;
3454 f -= 4;
3455
3456 ixt = ia + Hparam.xfirst - 1;
3457 iyt = ib + Hparam.yfirst - 1;
3458
3459 // xval1l = Hparam.xmin;
3460 // xval2l = Hparam.xmax;
3461 yval1l = Hparam.ymin;
3462 yval2l = Hparam.ymax;
3463
3464 xlab1l = gCurrentHist->GetXaxis()->GetXmin();
3465 xlab2l = gCurrentHist->GetXaxis()->GetXmax();
3466 if (Hoption.Logx) {
3467 if (xlab2l>0) {
3468 if (xlab1l>0) xlab1l = TMath::Log10(xlab1l);
3469 else xlab1l = TMath::Log10(0.001*xlab2l);
3470 xlab2l = TMath::Log10(xlab2l);
3471 }
3472 }
3473 ylab1l = gCurrentHist->GetYaxis()->GetXmin();
3474 ylab2l = gCurrentHist->GetYaxis()->GetXmax();
3475 if (Hoption.Logy) {
3476 if (ylab2l>0) {
3477 if (ylab1l>0) ylab1l = TMath::Log10(ylab1l);
3478 else ylab1l = TMath::Log10(0.001*ylab2l);
3479 ylab2l = TMath::Log10(ylab2l);
3480 }
3481 }
3482
3483 for (i = 1; i <= 4; ++i) {
3484 ixa = ixadd[i - 1];
3485 iya = iyadd[i - 1];
3486 Double_t xwid = gCurrentHist->GetXaxis()->GetBinWidth(ixt+ixa);
3487 Double_t ywid = gCurrentHist->GetYaxis()->GetBinWidth(iyt+iya);
3488
3489 // Compute the cell position in cartesian coordinates
3490 // and compute the LOG if necessary
3491 f[i*3 + 1] = gCurrentHist->GetXaxis()->GetBinLowEdge(ixt+ixa) + 0.5*xwid;
3492 f[i*3 + 2] = gCurrentHist->GetYaxis()->GetBinLowEdge(iyt+iya) + 0.5*ywid;
3493 if (Hoption.Logx) {
3494 if (f[i*3 + 1] > 0) f[i*3 + 1] = TMath::Log10(f[i*3 + 1]);
3495 else f[i*3 + 1] = Hparam.xmin;
3496 }
3497 if (Hoption.Logy) {
3498 if (f[i*3 + 2] > 0) f[i*3 + 2] = TMath::Log10(f[i*3 + 2]);
3499 else f[i*3 + 2] = Hparam.ymin;
3500 }
3501
3502 // Transform the cell position in the required coordinate system
3503 if (Hoption.System == kPOLAR) {
3504 f[i*3 + 1] = 360*(f[i*3 + 1] - xlab1l) / (xlab2l - xlab1l);
3505 f[i*3 + 2] = (f[i*3 + 2] - yval1l) / (yval2l - yval1l);
3506 } else if (Hoption.System == kCYLINDRICAL) {
3507 f[i*3 + 1] = 360*(f[i*3 + 1] - xlab1l) / (xlab2l - xlab1l);
3508 } else if (Hoption.System == kSPHERICAL) {
3509 f[i*3 + 1] = 360*(f[i*3 + 1] - xlab1l) / (xlab2l - xlab1l);
3510 f[i*3 + 2] = 360*(f[i*3 + 2] - ylab1l) / (ylab2l - ylab1l);
3511 } else if (Hoption.System == kRAPIDITY) {
3512 f[i*3 + 1] = 360*(f[i*3 + 1] - xlab1l) / (xlab2l - xlab1l);
3513 f[i*3 + 2] = (180 - dangle*2)*(f[i*3 + 2] - ylab1l) / (ylab2l - ylab1l) + dangle;
3514 }
3515
3516 // Get the content of the table. If the X index (ICX) is
3517 // greater than the X size of the table (NCX), that's mean
3518 // IGTABL tried to close the surface and in this case the
3519 // first channel should be used. */
3520 icx = ixt + ixa;
3521 if (icx > Hparam.xlast) icx = 1;
3522 f[i*3+3] = Hparam.factor*gCurrentHist->GetBinContent(icx, iyt + iya);
3523 if (Hoption.Logz) {
3524 if (f[i*3+3] > 0) f[i*3+3] = TMath::Log10(f[i*3+3]);
3525 else f[i*3+3] = Hparam.zmin;
3526 if (f[i*3+3] < Hparam.zmin) f[i*3+3] = Hparam.zmin;
3527 if (f[i*3+3] > Hparam.zmax) f[i*3+3] = Hparam.zmax;
3528 } else {
3529 f[i*3+3] = TMath::Max(Hparam.zmin, f[i*3+3]);
3530 f[i*3+3] = TMath::Min(Hparam.zmax, f[i*3+3]);
3531 }
3532
3533 // The colors on the surface can represent the content or the errors.
3534 // if (fSumw2.fN) t[i] = gCurrentHist->GetBinError(icx, iyt + iya);
3535 // else t[i] = f[i * 3 + 3];
3536 t[i] = f[i * 3 + 3];
3537 }
3538
3539 // Define the position of the colored contours for SURF3
3540 if (Hoption.Surf == 23) {
3541 for (i = 1; i <= 4; ++i) f[i * 3 + 3] = fRmax[2];
3542 }
3543
3544 if (Hoption.System == kCYLINDRICAL || Hoption.System == kSPHERICAL || Hoption.System == kRAPIDITY) {
3545 for (i = 1; i <= 4; ++i) {
3546 f[i*3 + 3] = (1 - rinrad)*((f[i*3 + 3] - Hparam.zmin) /
3547 (Hparam.zmax - Hparam.zmin)) + rinrad;
3548 }
3549 }
3550}
3551
3552////////////////////////////////////////////////////////////////////////////////
3553/// Draw surface in polar coordinates
3554///
3555/// \param[in] iordr order of variables (0 - R,PHI, 1 - PHI,R)
3556/// \param[in] na number of steps along 1st variable
3557/// \param[in] nb number of steps along 2nd variable
3558/// \param[in] chopt specific options
3559///
3560/// - `chopt` = 'BF' from BACK to FRONT
3561/// - `chopt` = 'FB' from FRONT to BACK
3562
3563void TPainter3dAlgorithms::SurfacePolar(Int_t iordr, Int_t na, Int_t nb, const char *chopt)
3564{
3565 /* Initialized data */
3566 static Int_t iface[4] = { 1,2,3,4 };
3567
3568 TView *view = gPad ? gPad->GetView() : nullptr;
3569 if (!view) {
3570 Error("SurfacePolar", "no TView in current pad");
3571 return;
3572 }
3573
3574 Int_t iphi, jphi, kphi, incr, nphi, iopt, iphi1, iphi2;
3575 Double_t f[12] /* was [3][4] */;
3576 Int_t i, j, incrr, ir1, ir2;
3577 Double_t z;
3578 Int_t ia, ib, ir, jr, nr, icodes[3]; // was icode[2]. One element more to differentiate front & back boxes from data
3579 Double_t tt[4];
3580 Double_t phi, ttt[4], xyz[12] /* was [3][4] */;
3581 ia = ib = 0;
3582
3583 if (iordr == 0) {
3584 jr = 1;
3585 jphi = 2;
3586 nr = na;
3587 nphi = nb;
3588 } else {
3589 jr = 2;
3590 jphi = 1;
3591 nr = nb;
3592 nphi = na;
3593 }
3594 if (fNaphi < nphi + 3) {
3595 fNaphi = nphi + 3;
3596 fAphi.resize(fNaphi);
3597 }
3598 if (fAphi.empty()) {
3599 Error("SurfacePolar", "failed to allocate array fAphi[%d]", fNaphi);
3600 fNaphi = 0;
3601 return;
3602 }
3603 iopt = 2;
3604 if (*chopt == 'B' || *chopt == 'b') iopt = 1;
3605
3606 // P R E P A R E P H I A R R A Y
3607 // F I N D C R I T I C A L S E C T O R S
3608 kphi = nphi;
3609 if (iordr == 0) ia = nr;
3610 if (iordr != 0) ib = nr;
3611 for (i = 1; i <= nphi; ++i) {
3612 if (iordr == 0) ib = i;
3613 if (iordr != 0) ia = i;
3614 (this->*fSurfaceFunction)(ia, ib, f, tt);
3615 if (i == 1) fAphi[0] = f[jphi - 1];
3616 fAphi[i - 1] = (fAphi[i - 1] + f[jphi - 1]) / (float)2.;
3617 fAphi[i] = f[jphi + 5];
3618 }
3619 view->FindPhiSectors(iopt, kphi, fAphi.data(), iphi1, iphi2);
3620
3621 // D R A W S U R F A C E
3622 icodes[2] = -1; // -1 for data, 0 for front a back boxes
3623 fEdgeIdx = 0; // constant since stacks are not (yet?) handled for surfaces
3624 incr = 1;
3625 iphi = iphi1;
3626L100:
3627 if (iphi > nphi) goto L300;
3628
3629 // F I N D O R D E R A L O N G R
3630 if (iordr == 0) {ia = nr; ib = iphi;}
3631 else {ia = iphi;ib = nr;}
3632
3633 (this->*fSurfaceFunction)(ia, ib, f, tt);
3634 phi = kRad*((f[jphi - 1] + f[jphi + 5]) / 2);
3635 view->FindNormal(TMath::Cos(phi), TMath::Sin(phi), 0, z);
3636 incrr = 1;
3637 ir1 = 1;
3638 if ((z <= 0 && iopt == 1) || (z > 0 && iopt == 2)) {
3639 incrr = -1;
3640 ir1 = nr;
3641 }
3642 ir2 = nr - ir1 + 1;
3643 // D R A W S U R F A C E F O R S E C T O R
3644 for (ir = ir1; incrr < 0 ? ir >= ir2 : ir <= ir2; ir += incrr) {
3645 if (iordr == 0) ia = ir;
3646 if (iordr != 0) ib = ir;
3647
3648 (this->*fSurfaceFunction)(ia, ib, f, tt);
3649 for (i = 1; i <= 4; ++i) {
3650 j = i;
3651 if (iordr != 0 && i == 2) j = 4;
3652 if (iordr != 0 && i == 4) j = 2;
3653 xyz[j*3 - 3] = f[jr + i*3 - 4]*TMath::Cos(f[jphi + i*3 - 4]*kRad);
3654 xyz[j*3 - 2] = f[jr + i*3 - 4]*TMath::Sin(f[jphi + i*3 - 4]*kRad);
3655 xyz[j*3 - 1] = f[i*3 - 1];
3656 ttt[j - 1] = tt[i - 1];
3657 }
3658 icodes[0] = ia;
3659 icodes[1] = ib;
3660 (this->*fDrawFace)(icodes, xyz, 4, iface, ttt);
3661 }
3662 // N E X T P H I
3663L300:
3664 iphi += incr;
3665 if (iphi == 0) iphi = kphi;
3666 if (iphi > kphi) iphi = 1;
3667 if (iphi != iphi2) goto L100;
3668 if (incr == 0) return;
3669 if (incr < 0) {
3670 incr = 0;
3671 goto L100;
3672 }
3673 incr = -1;
3674 iphi = iphi1;
3675 goto L300;
3676}
3677
3678////////////////////////////////////////////////////////////////////////////////
3679/// Draw surface in cylindrical coordinates
3680///
3681/// \param[in] iordr order of variables (0 - Z,PHI; 1 - PHI,Z)
3682/// \param[in] na number of steps along 1st variable
3683/// \param[in] nb number of steps along 2nd variable
3684/// \param[in] chopt specific options
3685///
3686/// - `chopt` = 'BF' from BACK to FRONT
3687/// - `chopt` = 'FB' from FRONT to BACK
3688
3689void TPainter3dAlgorithms::SurfaceCylindrical(Int_t iordr, Int_t na, Int_t nb, const char *chopt)
3690{
3691
3692
3693 /* Initialized data */