← Back to branch summary

~maddevelopers/mg5amcnlo/2.9.4

« back to all changes in this revision

Viewing changes to vendor/StdHEP/src/display/drawEvent.c

pass to v2.0.0

Show diffs side-by-side

added added

removed removed

Lines of Context:
vendor/StdHEP/src/display/drawEvent.c
 
1
/*******************************************************************************
 
2
*                                                                              *
 
3
* drawEvent.c -- Map the event - ie., the list of four-momentum vectors to an  *
 
4
*                array of points or segments for displaying on the spin widget *
 
5
*                                                                              *
 
6
* Copyright (c) 1991 Universities Research Association, Inc.                   *
 
7
* All rights reserved.                                                         *
 
8
*                                                                              *
 
9
* This material resulted from work developed under a Government Contract and   *
 
10
* is subject to the following license:  The Government retains a paid-up,      *
 
11
* nonexclusive, irrevocable worldwide license to reproduce, prepare derivative *
 
12
* works, perform publicly and display publicly by or for the Government,       *
 
13
* including the right to distribute to other Government contractors.  Neither  *
 
14
* the United States nor the United States Department of Energy, nor any of     *
 
15
* their employees, makes any warrenty, express or implied, or assumes any      *
 
16
* legal liability or responsibility for the accuracy, completeness, or         *
 
17
* usefulness of any information, apparatus, product, or process disclosed, or  *
 
18
* represents that its use would not infringe privately owned rights.           *
 
19
*                                                                              *
 
20
* Fermilab Nirvana GUI Library                                                 *
 
21
* June 17, 1991, August 22, 1991                                               *
 
22
*                                                                              *
 
23
* Written by Paul Lebrun & Mark Edel                                           *
 
24
*                                                                              *
 
25
*******************************************************************************/
 
26
static char SCCSID[] = "@(#)drawEvent.c 1.1     4/6/92";
 
27
#include <stdio.h>
 
28
#include <fcntl.h>
 
29
#include <stdlib.h>
 
30
#include <string.h>
 
31
#include <math.h>
 
32
#include <Xm/Xm.h>
 
33
#include "spin/Spin.h"
 
34
#include "util/stringUtils.h"
 
35
#include "phase.h"
 
36
#include "space.h"
 
37
#include "phaseP.h"
 
38
#include "drawEvent.h"
 
39
 
 
40
#define MAXPARTCOLOR 30         /* Size of particle color table */
 
41
#define EPSILON 0.1e-10
 
42
 
 
43
#define FULL_ON 0xffff
 
44
#define HALF_ON (0xffff/2)
 
45
#define OFF 0
 
46
 
 
47
#define NUMPINSPHI 8
 
48
#define NUMPINSTETA 7
 
49
#define MAXCHAR_HEPNAM 30  /* for use with hepnam_ */
 
50
 
 
51
extern float stdchg_(long *id);
 
52
 
 
53
static Pixel allocRGB(Screen *screen,
 
54
        unsigned short red, unsigned short green, unsigned short blue);
 
55
static Pixel particlePixel(int particleID);
 
56
static void detectorToSegment(SpaceWindow *window, SpinSegment *segdt, 
 
57
                                        SpinSegment *seg);
 
58
 
 
59
/*
 
60
** AllocateColors
 
61
**
 
62
** Generate the HEP color convention : This routines fills the 
 
63
** Global array ParticleColors, which is simply an list of color indices generated
 
64
** by calls to XAllocColor. 
 
65
**
 
66
** This convention is the following : 
 
67
**
 
68
**      The photon is always black.  If a color can not be allocated, Black 
 
69
**      is chosen.
 
70
**
 
71
**      For GrayScale Display, the charged leptons are white, the other 
 
72
**      neutrinos are black and the quarks and hadrons are gray. 
 
73
**
 
74
**      For Color display,  the charged leptons have the three fundamental 
 
75
**      RGB colors: electron is blue, muon is green and tau is red. 
 
76
**      The neutrinos are gray. The quarks are characterised a mixture
 
77
**      where one componenet of RGB is set to 0, the other two are 
 
78
**      half of the full intensity.  For instance, the up quark has 
 
79
**      no blue compoenet, and is half green, half red.
 
80
**      The Weak bosons and the gluons are made of a disgusting 
 
81
**      mixture of all RGB colors ( pink for the W andZ bosons, 
 
82
**      brown for the gluons).
 
83
**      The color of a hadron is set the by color of the heaviest quark 
 
84
**      in this hadron.
 
85
*/
 
86
void AllocateColors(Screen *screen)
 
87
{
 
88
    Display *display = DisplayOfScreen(screen);
 
89
    Pixel black = BlackPixelOfScreen(screen);
 
90
    Pixel white = WhitePixelOfScreen(screen);
 
91
    Pixel color;
 
92
    int i;
 
93
    XVisualInfo info;
 
94
 
 
95
    /* As a fallback, initialize all particle colors to black */
 
96
    for (i=0; i<MAXPARTCOLOR; i++)
 
97
        ParticleColors[i]=black;
 
98
 
 
99
    /* Assign colors based on type of screen available */
 
100
    if (XMatchVisualInfo(display, XScreenNumberOfScreen(screen),
 
101
            DefaultDepthOfScreen(screen), PseudoColor, &info)) {
 
102
        /*
 
103
        ** Screen can handle PseudoColor (color lookup table)
 
104
        */
 
105
        /* The Up quark */
 
106
        ParticleColors[1] = allocRGB(screen, FULL_ON, FULL_ON, OFF);
 
107
 
 
108
        /* The Down quark */
 
109
        ParticleColors[2] = allocRGB(screen, FULL_ON, HALF_ON, OFF);
 
110
 
 
111
        /* The Strange quark */
 
112
        ParticleColors[3] = allocRGB(screen, OFF, FULL_ON, FULL_ON);
 
113
 
 
114
        /* The Charm quark */
 
115
        ParticleColors[4] = allocRGB(screen, OFF, FULL_ON, HALF_ON);
 
116
 
 
117
        /* The B  quark */
 
118
        ParticleColors[5] = allocRGB(screen, FULL_ON, OFF, FULL_ON);
 
119
 
 
120
        /* The Top quark */
 
121
        ParticleColors[6] = allocRGB(screen, HALF_ON, OFF, FULL_ON);
 
122
 
 
123
        /* The electron  */
 
124
        ParticleColors[11] = allocRGB(screen, OFF, OFF, FULL_ON);
 
125
 
 
126
        /* The muon */
 
127
        ParticleColors[13] = allocRGB(screen, OFF, FULL_ON, OFF);
 
128
 
 
129
        /* The Tau  */
 
130
        ParticleColors[15] = allocRGB(screen, FULL_ON, OFF, OFF);
 
131
 
 
132
        /* Neutrinos  */
 
133
        color = allocRGB(screen, FULL_ON/3, FULL_ON/3, FULL_ON/3);
 
134
        ParticleColors[16] = ParticleColors[14] = ParticleColors[12] = color;
 
135
 
 
136
        /* The photon quark */
 
137
        ParticleColors[22] = black;
 
138
 
 
139
        /* The W boson */
 
140
        ParticleColors[24] = allocRGB(screen, FULL_ON, HALF_ON, HALF_ON);
 
141
 
 
142
        /* The Z Boson */
 
143
        ParticleColors[23] = allocRGB(screen, FULL_ON, FULL_ON, HALF_ON);
 
144
 
 
145
        /* Gluons */
 
146
        color = allocRGB(screen, FULL_ON, HALF_ON, 2*FULL_ON/3);
 
147
        ParticleColors[21] = ParticleColors[9] = color;
 
148
 
 
149
    } else if (XMatchVisualInfo(display, XScreenNumberOfScreen(screen),
 
150
            DefaultDepthOfScreen(screen), GrayScale, &info)) {
 
151
        /*
 
152
        ** Screen can't handle color but can handle GrayScale
 
153
        */
 
154
        /* Quarks */
 
155
        color = allocRGB(screen, OFF, FULL_ON, FULL_ON);
 
156
        for (i=1; i<7; i++)
 
157
            ParticleColors[i] = color;
 
158
        ParticleColors[9] = color;
 
159
        ParticleColors[21]= color;
 
160
 
 
161
        /* Leptons */
 
162
        ParticleColors[11] = white;
 
163
        ParticleColors[13] = white;
 
164
        ParticleColors[15] = white;
 
165
    }
 
166
}
 
167
/*
 
168
** Produce a bunch of Radial coordinates to build a small icon describing the 
 
169
** this vertex.
 
170
*/
 
171
void SetPinVertex()
 
172
{
 
173
  int i,j;
 
174
  double pi, dphi, dteta, phi, teta;
 
175
  
 
176
  pi = acos(-1.0);
 
177
  dphi = 2.0 * pi / NUMPINSPHI; dteta = pi / NUMPINSTETA;
 
178
  for (i=0, phi = 0. ; i < NUMPINSTETA; i++, phi += dphi) {
 
179
    for (j=0, teta = 0. ; j < NUMPINSPHI; j++, teta += dteta) {
 
180
      XPinVertexPoints[j][i] = sin(teta) * cos(phi);
 
181
      YPinVertexPoints[j][i] = sin(teta) * sin(phi);
 
182
      ZPinVertexPoints[j][i] = cos(teta);
 
183
    }
 
184
  }
 
185
}
 
186
 
 
187
/*
 
188
** DrawEvent
 
189
**
 
190
** Map an event to 3D segments for the spin widget.  From a 4-momentum vector,
 
191
** (i.e., the 3-momentum vector and the unique identifier from which a physical
 
192
** mass can be obtained), An array of 3-dimensional vector and track type is
 
193
** deduced).  The structure PhaseEvent is described in phase.h.  The setting
 
194
** of the displayMode toggle on the event display panel determines how the
 
195
** event is mapped onto the screen:
 
196
** 
 
197
** BY_MOMENTUM:         The length is proportional the  |P|, the norm of the 3-d
 
198
**                      momentum vector. This is the simplest and default option.
 
199
**
 
200
** BY_PT:               The length of the vector is proportional to Pt, the 
 
201
**                      transverse momentum. This trick greatly enhances 
 
202
**                      the "high Pt" Physics in the event.
 
203
**
 
204
** BY_RAPIDITY:         The length of the vector is propotional to the rapidity 
 
205
**                      of the particle, as defined in the PDG bible.
 
206
**
 
207
** BY_PSEUDORAPIDITY:   The length of the vector is proportional to the pseudo
 
208
**                      rapidity, i.e., the mass of the particle is 
 
209
**                      neglected.
 
210
**
 
211
** BY_USER_ROUTINE :    The user function specified in the DisplayEvent call
 
212
**                      is called to determine what the new vector ought to be... 
 
213
*/
 
214
void DrawEvent(StdHepWindow *window, int setScale)
 
215
{
 
216
    int i, maxIndex, nSegments, nns, npvis;
 
217
    int doPoints;
 
218
    double scale, innerScale, length, maxLength, outVec[3];
 
219
    SpinSegment *segment, *segments, *segdt, temp;
 
220
    SpinPoint *point, *points;
 
221
    PhaseParticle *p, *particles = window->event.particles;
 
222
    int nStable = 0, nParticles = window->event.nParticles;
 
223
    SpaceVertex *pinv, *v;
 
224
    char nameText[MAXCHAR_HEPNAM], statText[255];
 
225
    XmString cStatText;
 
226
    /* There is a very wierd Ultrix compiler problem concerning this statement:*/
 
227
    /*    PhaseEvent *event; */
 
228
    static char *unitText[6] = {"p", "Pt", "Eta", "Eta", "value", "Eta-Pt"};
 
229
    XmFontList fontList;
 
230
    PhaseWindow *winp;
 
231
    SpaceWindow *wins;
 
232
    ParaWindow *winpa;
 
233
    
 
234
    /* If a Space window, find out the number of vertices in this event, 
 
235
       so we can allocate more space to draw short segments to represent 
 
236
       a vertex.. */
 
237
       
 
238
    if (window->type == STDHEP_SPACE) {
 
239
       wins = (SpaceWindow *) window;
 
240
    /* allocate space for the maximum possible number of segments or points */
 
241
       nns = NUMPINSTETA * NUMPINSPHI * wins->numRealVertices + nParticles;
 
242
       nns += wins->nDetectorSegments;
 
243
       segments = (SpinSegment *)XtMalloc(sizeof(SpinSegment) * nns);
 
244
       doPoints = False;
 
245
     } else if (window->type == STDHEP_PARA) {
 
246
        winpa = (ParaWindow *) window;
 
247
        segments = (SpinSegment *)XtMalloc(sizeof(SpinSegment) * nParticles);
 
248
        doPoints = False;
 
249
     } else {
 
250
        winp =(PhaseWindow *) window;
 
251
        doPoints = (winp->vectorLength == 0.);
 
252
        if (doPoints)
 
253
            points = (SpinPoint *)XtMalloc(sizeof(SpinPoint) * nParticles);
 
254
        else
 
255
           segments = (SpinSegment *)XtMalloc(sizeof(SpinSegment) * nParticles);
 
256
     }
 
257
    
 
258
    /* Loop over all of the particles, translating them into 3-D line segments
 
259
       or points for the Spin widget.  Also find the particle of maximum value
 
260
       (for setting the scale and printing statistics), and the number of stable
 
261
       particles (for printing statistics later) */
 
262
    nSegments = 0;
 
263
    npvis = 0;
 
264
    maxLength = -1.;
 
265
    if (doPoints) {
 
266
        point = points;
 
267
        for (i=0, p=particles; i<nParticles; i++, p++) {
 
268
            if (p->stable)
 
269
                nStable++;
 
270
            if (ParticleToSegment(winp, p, &temp, &length)) {
 
271
                point->x = temp.x1;
 
272
                point->y = temp.y1;
 
273
                point->z = temp.z1;
 
274
                point->pixel = temp.pixel;
 
275
                point++;
 
276
                nSegments++;
 
277
                npvis++;
 
278
                if (length > maxLength) {
 
279
                    maxLength = length;
 
280
                    maxIndex = i;
 
281
                }
 
282
            }
 
283
        }
 
284
        /* Display the points */
 
285
        SpinSetPoints(window->spin, points, nSegments);
 
286
        XtFree((char *) points);
 
287
    } else {
 
288
      if (window->type == STDHEP_PHASE) {
 
289
        segment=segments;
 
290
        for (i=0, p=particles; i<nParticles; i++, p++) {
 
291
            if (p->stable)
 
292
                nStable++;
 
293
            if (ParticleToSegment(winp, p, segment, &length)) {
 
294
                segment++;
 
295
                nSegments++;
 
296
                npvis++;
 
297
                if (length > maxLength) {
 
298
                    maxLength = length;
 
299
                    maxIndex = i;
 
300
                }
 
301
            }
 
302
        }
 
303
       } else if (window->type == STDHEP_PARA) {
 
304
         segment=segments;
 
305
         for (i=0, p=particles; i<nParticles; i++, p++) {
 
306
            if (p->stable)
 
307
                nStable++;
 
308
            if (ParticleToParaSegment(winpa, p, segment, &length)) {
 
309
                segment++;
 
310
                nSegments++;
 
311
                npvis++;
 
312
                if (length > maxLength) {
 
313
                    maxLength = length;
 
314
                    maxIndex = i;
 
315
                }
 
316
            }
 
317
          }
 
318
       } else {            
 
319
           segment=segments;
 
320
           maxLength = 0.;
 
321
           for (i=0, p=particles, v=wins->vertices; 
 
322
                             i<nParticles; i++, p++, v++) {
 
323
            if (p->stable)
 
324
                nStable++;
 
325
            if (TrackToSegment(wins, p, v, segment, &length)) {
 
326
                segment++;
 
327
                nSegments++;
 
328
                npvis++;
 
329
                if (length > maxLength) {
 
330
                    maxLength = length;
 
331
                    maxIndex = i;
 
332
                }
 
333
            }
 
334
           }
 
335
           if (wins->nDetectorSegments > 0) {
 
336
            for (i=0, segdt = wins->detectorSegments;
 
337
                  i < wins->nDetectorSegments;
 
338
                   i++, segdt++, segment++, nSegments++)  
 
339
                  detectorToSegment(wins, segdt, segment);
 
340
           } 
 
341
           if (wins->showVertices) {
 
342
           for (i=0, pinv=wins->realVertices;
 
343
                 i < wins->numRealVertices; i++, pinv++){
 
344
                  if (fabs((double) pinv->z) < wins->currentScale) {
 
345
                    VertexToSegment(wins, pinv, &segment);
 
346
                     nSegments += NUMPINSPHI *NUMPINSTETA;
 
347
                  } 
 
348
             }
 
349
           }
 
350
         }
 
351
        }    
 
352
        /* Display the line segments with the spin widget */
 
353
        SpinSetSegments(window->spin, segments, nSegments);
 
354
        XtFree((char *) segments);
 
355
        /* 
 
356
        * Currently, the spin widget Perspective resource is buggy, 
 
357
        * turn if off.. 
 
358
        */
 
359
        SET_ONE_RSRC(window->spin, XmNperspectiveOn, FALSE);
 
360
    
 
361
    /* If requested, set the scale to show the data most effectively */
 
362
    if (setScale) {
 
363
        SpinSetScale(window->spin, (maxLength > 0.) ? 1./maxLength : 1.);
 
364
        DrawScale(window);
 
365
    }
 
366
    
 
367
    /* Fill in the statistics window */
 
368
    if (npvis == 0) {
 
369
        sprintf(statText, "%d tracks, %d are stable.\nNo Tracks are visible",
 
370
                nParticles, nStable);
 
371
    } else {
 
372
        if ((window->type == STDHEP_PHASE) || (window->type == STDHEP_PARA)) { 
 
373
            GET_ONE_RSRC(window->statsLabel, XmNfontList, &fontList);
 
374
            if (CharsetAvailable("greek", fontList)) {
 
375
                hepnmg_(&particles[maxIndex].id, nameText);
 
376
            } else {
 
377
                hepnam_(&particles[maxIndex].id, nameText, MAXCHAR_HEPNAM);
 
378
                 *strchr(nameText, ' ') = '\0';
 
379
            }
 
380
            if (window->type == STDHEP_PHASE) {
 
381
            sprintf(statText,
 
382
"%d tracks, %d are stable,\n%d are now displayed\nTrack of max. %s is a(n) %s",
 
383
                nParticles, nStable, nSegments,
 
384
                unitText[winp->displayMode], nameText);
 
385
            } else { 
 
386
            sprintf(statText,
 
387
"%d tracks, %d are stable,\n%d are now displayed\nTrack of max. %s is a(n) %s",
 
388
                nParticles, nStable, nSegments,
 
389
                unitText[5], nameText);
 
390
            }
 
391
          } else {
 
392
            if (wins->showVertices)
 
393
              nns = nSegments
 
394
                 - (NUMPINSPHI * NUMPINSTETA * wins->numRealVertices);
 
395
            else nns = nSegments; 
 
396
            sprintf(statText,
 
397
"%d tracks, %d are stable,\n%d are now displayed\n %d distinct vertices found",
 
398
                nParticles, nStable, nns, wins->numRealVertices);
 
399
          }
 
400
    }
 
401
    cStatText = MultiFontString(statText);
 
402
    SET_ONE_RSRC(window->statsLabel, XmNlabelString, cStatText);
 
403
    XmStringFree(cStatText);
 
404
}
 
405
 
 
406
int ParticleToSegment(PhaseWindow *window, PhaseParticle *p,
 
407
                      SpinSegment *seg, double *length)
 
408
{
 
409
    int outColor;
 
410
    double scale, innerScale, outVec[3];
 
411
 
 
412
    if (!ParticleVisible((StdHepWindow *) window, p->id, p->stable))
 
413
        return False;
 
414
 
 
415
    switch (window->displayMode) {
 
416
      case BY_MOMENTUM:
 
417
        *length = ParticleMomentum(p->px, p->py, p->pz);
 
418
        break;
 
419
      case BY_PT:
 
420
        *length = ParticlePT(p->px, p->py);
 
421
        break;
 
422
      case BY_RAPIDITY:
 
423
        *length = ParticleRapidity(p->px, p->py, p->pz, p->mass);
 
424
        break;
 
425
      case BY_PSEUDORAPIDITY:
 
426
        *length = ParticlePseudorapidity(p->px, p->py, p->pz);
 
427
        break;
 
428
/*       case BY_USER_ROUTINE:
 
429
        (*window->userMapProc)(window->event, outVec, &outColor);
 
430
        *length = sqrt(pow(outVec[0],2) + pow(outVec[1],2) + pow(outVec[2],2));
 
431
        break; */
 
432
    }
 
433
    if (window->displayMode == BY_USER_ROUTINE) {
 
434
        seg->x1 = outVec[0];
 
435
        seg->y1 = outVec[1];
 
436
        seg->z1 = outVec[2];
 
437
        seg->x2 = outVec[0] * (1. - window->vectorLength);
 
438
        seg->y2 = outVec[0] * (1. - window->vectorLength);
 
439
        seg->z2 = outVec[0] * (1. - window->vectorLength);
 
440
        seg->pixel = ParticleColors[outColor];
 
441
    } else {
 
442
        scale = *length / ParticleMomentum(p->px, p->py, p->pz);
 
443
        innerScale = scale * (1. - window->vectorLength);
 
444
        seg->x1 = p->px * scale;
 
445
        seg->y1 = p->py * scale;
 
446
        seg->z1 = p->pz * scale;
 
447
        seg->x2 = p->px * innerScale;
 
448
        seg->y2 = p->py * innerScale;
 
449
        seg->z2 = p->pz * innerScale;
 
450
        seg->pixel = particlePixel(p->id);
 
451
    }
 
452
    return True;
 
453
}
 
454
 
 
455
int ParticleToParaSegment(ParaWindow *window, PhaseParticle *p,
 
456
                      SpinSegment *seg, double *length)
 
457
{
 
458
    if (!ParticleVisible((StdHepWindow *) window, p->id, p->stable))
 
459
        return False;
 
460
 
 
461
     seg->x1 = 0.;
 
462
     seg->y1 = 0.;
 
463
     seg->x2 = p->px;
 
464
     seg->y2 = p->py;
 
465
     if (p->pz < 0.) 
 
466
         seg->z1 = window->currentTranslz - 
 
467
               ( ParticleRapidity(p->px, p->py, p->pz, p->mass) * 
 
468
                window->currentRapToPt);
 
469
     else 
 
470
         seg->z1 = window->currentTranslz + 
 
471
               ( ParticleRapidity(p->px, p->py, p->pz, p->mass) * 
 
472
                window->currentRapToPt);
 
473
     seg->z2 = seg->z1;
 
474
     seg->pixel = particlePixel(p->id);
 
475
     *length = ParticlePT(p->px, p->py);
 
476
    return True;
 
477
}
 
478
 
 
479
int TrackToSegment(SpaceWindow *window,  PhaseParticle *p, SpaceVertex *v,
 
480
                      SpinSegment *seg, double *length)
 
481
{
 
482
        double scaleZ, scaleT;
 
483
        
 
484
        if (!ParticleVisible((StdHepWindow *)window, p->id, p->stable))
 
485
        return False;
 
486
        
 
487
        scaleZ = window->currentMomToSpace;
 
488
        scaleT = window->currentMomToSpace * window->currentLongToTr;
 
489
        seg->x1 = window->currentLongToTr * v->x;
 
490
        seg->y1 = window->currentLongToTr * v->y;
 
491
        seg->z1 = v->z;
 
492
        /*
 
493
        ** Check now that we are in the visible range
 
494
        **
 
495
        */
 
496
        if (fabs(seg->z1) > window->currentScale) return False;
 
497
        /*
 
498
        ** Now translate 
 
499
        */
 
500
        seg->x1 = seg->x1 + window->currentLongToTr * window->currentTransl[0];
 
501
        seg->y1 = seg->y1 + window->currentLongToTr * window->currentTransl[1];
 
502
        seg->z1 = seg->z1 + window->currentTransl[2];
 
503
        /* 
 
504
        ** The end of the tracks is set by the momentum
 
505
        */
 
506
        seg->x2 = seg->x1 + scaleT * p->px;
 
507
        seg->y2 = seg->y1 + scaleT * p->py;
 
508
        seg->z2 = seg->z1 + scaleZ * p->pz;
 
509
        *length = sqrt ((seg->x2 * seg->x2) + (seg->y2 * seg->y2) + 
 
510
                         (seg->z2 * seg->z2));
 
511
        seg->pixel = particlePixel(p->id);
 
512
        return True;
 
513
}
 
514
 
 
515
void VertexToSegment(SpaceWindow *window,  SpaceVertex *v,
 
516
                      SpinSegment **seg)
 
517
{
 
518
        double scaleZ, scaleT;
 
519
        int i,j;
 
520
        SpinSegment *segt;
 
521
        Pixel black;
 
522
        Screen *screen;
 
523
        
 
524
        screen = XtScreen(window->shell);
 
525
        black = BlackPixelOfScreen(screen);      
 
526
        scaleZ = PINSIZE * window->maxTransl[2];
 
527
/*      scaleT = PINSIZE * window->maxTransl[0]; */
 
528
        segt = *seg;
 
529
        for (i=0; i < NUMPINSTETA; i++) {
 
530
          for (j=0; j < NUMPINSPHI; j++) {
 
531
            segt->x1 = window->currentLongToTr * v->x;
 
532
            segt->y1 = window->currentLongToTr * v->y;
 
533
            segt->z1 = v->z;
 
534
        /*
 
535
        ** Now translate 
 
536
        */
 
537
            segt->x1 = segt->x1 + 
 
538
               window->currentLongToTr * window->currentTransl[0];
 
539
            segt->y1 = segt->y1 + 
 
540
               window->currentLongToTr * window->currentTransl[1];
 
541
            segt->z1 = segt->z1 + 
 
542
               window->currentTransl[2];
 
543
            /*
 
544
            ** The end of the pins.. 
 
545
            */
 
546
            segt->x2 = segt->x1 + scaleZ * XPinVertexPoints[j][i];
 
547
            segt->y2 = segt->y1 + scaleZ * YPinVertexPoints[j][i];
 
548
            segt->z2 = segt->z1 + scaleZ * ZPinVertexPoints[j][i];
 
549
            segt->pixel = black;
 
550
            segt++;
 
551
            }
 
552
         }
 
553
         *seg = segt;
 
554
}
 
555
 
 
556
static void detectorToSegment(SpaceWindow *window, SpinSegment *segdt, 
 
557
                                        SpinSegment *seg)
 
558
{
 
559
        int i,j;
 
560
        Pixel black;
 
561
        Screen *screen;
 
562
        
 
563
        screen = XtScreen(window->shell);
 
564
        black = BlackPixelOfScreen(screen);      
 
565
        seg->pixel = black;
 
566
        seg->x1 = window->currentLongToTr * segdt->x1;
 
567
        seg->y1 = window->currentLongToTr * segdt->y1;
 
568
        seg->z1 = segdt->z1;
 
569
        seg->x2 = window->currentLongToTr * segdt->x2;
 
570
        seg->y2 = window->currentLongToTr * segdt->y2;
 
571
        seg->z2 = segdt->z2;
 
572
        /*
 
573
        ** Now translate 
 
574
        */
 
575
        seg->x1 = seg->x1 + 
 
576
               window->currentLongToTr * window->currentTransl[0];
 
577
        seg->y1 = seg->y1 + 
 
578
               window->currentLongToTr * window->currentTransl[1];
 
579
        seg->z1 = seg->z1 + 
 
580
               window->currentTransl[2];
 
581
        seg->x2 = seg->x2 + 
 
582
               window->currentLongToTr * window->currentTransl[0];
 
583
        seg->y2 = seg->y2 + 
 
584
               window->currentLongToTr * window->currentTransl[1];
 
585
        seg->z2 = seg->z2 + 
 
586
               window->currentTransl[2];
 
587
        return;
 
588
}
 
589
 
 
590
 
 
591
        
 
592
void SetDisplayMode(PhaseWindow *window, int newMode)
 
593
{
 
594
    int changeScale, oldMode = window->displayMode;
 
595
    
 
596
    changeScale = (((newMode == BY_MOMENTUM || newMode == BY_PT) &&
 
597
          (oldMode != BY_MOMENTUM && oldMode != BY_PT)    ) ||
 
598
         ((newMode == BY_RAPIDITY || newMode == BY_PSEUDORAPIDITY) &&
 
599
          (oldMode != BY_RAPIDITY && oldMode != BY_PSEUDORAPIDITY)    ) ||
 
600
         (newMode == BY_USER_ROUTINE && oldMode != BY_USER_ROUTINE));
 
601
    window->displayMode = newMode;
 
602
    DrawEvent((StdHepWindow *) window, changeScale);
 
603
}
 
604
 
 
605
void DrawScale(StdHepWindow *window)
 
606
{
 
607
    Widget scaleWidget = window->scaleArea;
 
608
    Window scaleWindow = XtWindow(scaleWidget);
 
609
    Display *display = XtDisplay(scaleWidget);
 
610
    double scale, scaleAreaWidth, scaleBarLength = 1., barMultiplier = 2.;
 
611
    double tempMom;
 
612
    short widgetWidth, spinWidth, widgetHeight;
 
613
    int leftX, rightX, topY, midY, bottomY;
 
614
    GC gc = window->scaleGC;
 
615
    char scaleText[100];
 
616
    XmString cScaleText;
 
617
    SpaceWindow *wins;
 
618
    double scaleVertex;
 
619
    PhaseWindow *winp;
 
620
    ParaWindow *winpa;
 
621
 
 
622
    /* Get the scale factor for translating spin coordinates to pixels */
 
623
    GET_ONE_RSRC(window->spin, XmNwidth, &spinWidth);
 
624
    scale = (1./(SpinGetScale(window->spin)/2.)) / (double)spinWidth;
 
625
    
 
626
    /* Calculate the width of the scale draw area in spin coordinates */
 
627
    GET_ONE_RSRC(scaleWidget, XmNwidth, &widgetWidth);
 
628
    scaleAreaWidth = widgetWidth * scale;
 
629
    
 
630
    /* choose a sensible length for the scale bar, such as 1, 5, 10, 50 etc. */
 
631
    if (scaleAreaWidth < 1.) {
 
632
        while (scaleBarLength > scaleAreaWidth) {
 
633
            scaleBarLength /= barMultiplier;
 
634
            barMultiplier = (barMultiplier == 2.) ? 5. : 2.;
 
635
        }
 
636
    } else {
 
637
        while (scaleBarLength < scaleAreaWidth) {
 
638
            barMultiplier = (barMultiplier == 2.) ? 5. : 2.;
 
639
            scaleBarLength *= barMultiplier;
 
640
        }
 
641
        scaleBarLength /= barMultiplier;
 
642
    }
 
643
    
 
644
    /* draw the scale bar in the middle of the scale area widget */
 
645
    GET_ONE_RSRC(scaleWidget, XmNheight, &widgetHeight);
 
646
    leftX = ((scaleAreaWidth - scaleBarLength) / 2.) / scale;
 
647
    rightX = leftX + scaleBarLength / scale;
 
648
    topY = widgetHeight/4;
 
649
    midY = widgetHeight/2;
 
650
    bottomY = midY + widgetHeight/4;
 
651
    XClearWindow(display, scaleWindow);
 
652
    XDrawLine(display, scaleWindow, gc, leftX, topY, leftX, bottomY);
 
653
    XDrawLine(display, scaleWindow, gc, rightX, topY, rightX, bottomY);
 
654
    XDrawLine(display, scaleWindow, gc, leftX, midY, rightX, midY);
 
655
    
 
656
    /* put in the text to label the scale */
 
657
    if (window->type == STDHEP_PHASE) {
 
658
      winp = (PhaseWindow *) window;
 
659
      if (winp->displayMode == BY_MOMENTUM || winp->displayMode == BY_PT)
 
660
        sprintf(scaleText, "%G Gev/c", scaleBarLength);
 
661
      else
 
662
        sprintf(scaleText, "%G Unit Eta", scaleBarLength);
 
663
    } else if (window->type == STDHEP_PARA) {
 
664
      winpa = (ParaWindow *) window;
 
665
      scaleVertex = scaleBarLength / winpa->currentRapToPt;
 
666
      sprintf (scaleText, "%G eta = %G Gev/c",scaleBarLength, scaleVertex);
 
667
    
 
668
    } else {
 
669
      wins = (SpaceWindow *) window;
 
670
      tempMom = wins->currentMomToSpace;
 
671
      if( tempMom == 0.0) tempMom = EPSILON;
 
672
      scaleVertex = scaleBarLength / tempMom;
 
673
      sprintf (scaleText, "%G cm = %G Gev/c",scaleBarLength, scaleVertex);
 
674
    } 
 
675
    cScaleText = MKSTRING(scaleText);
 
676
    SET_ONE_RSRC(window->scaleLabel, XmNlabelString, cScaleText);
 
677
    XmStringFree(cScaleText);
 
678
}
 
679
 
 
680
 
 
681
static Pixel allocRGB(Screen *screen,
 
682
        unsigned short red, unsigned short green, unsigned short blue)
 
683
{
 
684
    Display *display = DisplayOfScreen(screen);
 
685
    XColor color;
 
686
    
 
687
    color.red = red;
 
688
    color.green = green;
 
689
    color.blue = blue;
 
690
    if (XAllocColor(display, DefaultColormapOfScreen(screen), &color))
 
691
        return color.pixel;
 
692
    else 
 
693
        return BlackPixelOfScreen(screen);
 
694
}
 
695
 
 
696
int ParticleVisible(StdHepWindow *window, long id, int stable)
 
697
{
 
698
    long absID = abs(id);
 
699
    float charge = stdchg_(&id);
 
700
    
 
701
    if (stable && !window->stable)
 
702
        return False;
 
703
    if (!stable && !window->unstable)
 
704
        return False;
 
705
        
 
706
    if ((fabs(charge) < 0.1) && !window->neutral)
 
707
        return False;
 
708
    if ((fabs(charge) > 0.1) && !window->charged)
 
709
        return False;
 
710
     
 
711
    if (absID == 11 && !window->electrons)
 
712
        return False;
 
713
    if (absID == 13 && !window->muons)
 
714
        return False;
 
715
    if ((absID == 12 || absID == 14 || absID == 16) && !window->neutrinos)
 
716
        return False;
 
717
    if ((absID <=  6 || absID == 9 || absID == 21) && !window->quarks)
 
718
        return False;
 
719
    if ((absID == 23 || absID == 24) && !window->wz)
 
720
        return False; 
 
721
    if (absID == 22 && !window->gammas)
 
722
        return False;
 
723
    if (absID >= 100 && !window->hadrons)
 
724
        return False;
 
725
    return True;
 
726
}
 
727
 
 
728
static Pixel particlePixel(int particleID)
 
729
{
 
730
    return ParticleColors[ParticleColorIndex(particleID)];
 
731
}
 
732
 
 
733
int ParticleColorIndex(int particleID)
 
734
{
 
735
    long absID = abs(particleID);
 
736
    long i, aid, kq[3], kqj;
 
737
    double akqx, akq3, akq2, akq1, akqj; 
 
738
    
 
739
    /*
 
740
    ** For non hadrons, return the particle id as the color
 
741
    */
 
742
    if (absID < 100) {
 
743
        if (absID > (MAXPARTCOLOR - 1))
 
744
            return MAXPARTCOLOR - 1;
 
745
        return absID;
 
746
    }
 
747
    
 
748
    /*
 
749
    ** For hadrons, return the largest digit? (ask paul...)
 
750
    */
 
751
    akqj = ((double)absID)/10.;
 
752
 
 
753
    aid = (long)akqj;
 
754
    akq1 = akqj/10.;
 
755
    kq[0] = aid - ((long)akq1) * 10;
 
756
 
 
757
    aid = (long)akq1;
 
758
    akq2 = akq1/10.;
 
759
    kq[1] = aid - ((long)akq2) * 10;
 
760
 
 
761
    aid = (long)akq2;
 
762
    akq3 = akq2/10.;
 
763
    kq[2] = aid - ((long)akq3) * 10;
 
764
 
 
765
    kqj = 0;
 
766
    for (i=0; i<3; i++)
 
767
        if (kq[i] > kqj) kqj = kq[i];  
 
768
 
 
769
    if (kqj > (MAXPARTCOLOR-1))
 
770
        return MAXPARTCOLOR-1;
 
771
    else
 
772
        return kqj;
 
773
}
 
774
 
 
775
double ParticleMomentum(double px, double py, double pz)
 
776
{
 
777
    double pTmp = sqrt(SQR(px) + SQR(py) + SQR(pz));
 
778
    if (pTmp == 0.0) return EPSILON;
 
779
    else return pTmp;
 
780
}
 
781
 
 
782
double ParticlePT(double px, double py)
 
783
{
 
784
    return sqrt(SQR(px) + SQR(py));
 
785
}
 
786
 
 
787
double ParticleRapidity(double px, double py, double pz, double mass)
 
788
{
 
789
    double p, pt, e;
 
790
    
 
791
    p = ParticleMomentum(px, py, pz);
 
792
    e = sqrt(p*p + mass*mass);
 
793
    return fabs(0.5*log(fabs((e+pz+0.0001)/(e-pz+0.0000001))));
 
794
}
 
795
 
 
796
double ParticlePseudorapidity(double px, double py, double pz)
 
797
{
 
798
    double p, pt, teta;
 
799
    
 
800
    p = ParticleMomentum(px, py, pz);
 
801
    pt = ParticlePT(px,py);
 
802
    teta = acos(pz/(p+ 0.00001));
 
803
    return fabs(log(fabs(tan(0.5*teta))));
 
804
}