-
Notifications
You must be signed in to change notification settings - Fork 4
/
hddsCommon.cpp
2305 lines (2159 loc) · 66.9 KB
/
hddsCommon.cpp
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
/* HDDS Common Classes
*
* Author: [email protected]
*
* Original version - Richard Jones, January 6, 2006.
* Revision 1.1 - Richard Jones, September 10, 2013.
*
* Revision 1.1 Notes:
* -------------------
* 1. Updated to HDDS-1.1 which adds support for geometry layers, using
* the geometry_layer="int" attribute of the object placement tags.
*
* Original Notes:
* ---------------
* 1. The HDDS specification is an xml document in the standard W3C
* schema namespace http://www.gluex.org/hdds, as described by the
* HDDS-1_0.xsd schema document.
* 2. Access to the xml source is through the industry-standard DOM interface.
* 3. The code has been tested with the xerces-c DOM implementation from
* Apache, and is intended to be used with the xerces-c library.
* 4. The common classes were created originally as tools for the hdds-geant
* converter to write Fortran geometry definition code for the Geant3
* simulation. Later they were reused for other simulation packages,
* but occasionally the original purpose is still visible in comments
* and some of the implementation details.
*
* Implementation:
* ---------------
* Most of the translation was straight-forward, but there were a couple
* of tricky cases where decisions had to be made. I think that these
* choices should work out in most cases. If not, further tuning of the
* algorithm will be necessary. Here are the tricky cases.
*
* 1. When to use divisions instead of placing multiple copies.
*
* Most of the time when a <mpos...> command appears it can be translated
* into a division of the mother volume in Geant. This is good to do
* because it makes both more compact description and is more efficient
* at tracking time. The difficulty here is that there is no easy way
* to check if the contents fit entirely inside the division. This is
* not a problem in the HDDS geometry description because the <mpos...>
* command is only for positioning, and makes no statement about what
* slice of the mother it occupies. But it is a problem for Geant because
* contents of divisions have to fit inside the division. This can
* happen at any time, but it most frequently happens when the object
* is being rotated before placement, as in the case of stereo layers
* in a drift chamber. My solution is to make a strict set of rules
* that are required before hdds-geant will create a division in response
* to a <mpos...> command, and do individual placement by default.
* The rules for creation of divisions are as follows:
* (a) the <composition> command must have a solid container, either
* via the envelope="..." attribute or by itself being a division.
* (b) the <mpos...> command must be alone inside its <composition>
* (c) the kind of shape of the container must be compatible with the
* <mpos...> type, eg. <mposPhi> would work if its container is
* a "tubs" (or division theroef) but not if it is a "box".
* (d) for <mposPhi> the impliedRot attribute must be "true".
* (e) the rot="..." attribute must be zeros or missing.
* The last condition is not logically necessary for it to work, but it
* avoids the problems that often occur with rotated placements failing
* to fit inside the division. To bypass this limitation and place
* rotated volumes inside divisions, one can simply create a new volume
* as a <composition> into which the content is placed with rotation,
* and then place the new volume with the <mpos...> command without rot.
*
* 2. How to recognize which media contain magnetic fields.
*
* There is was originally no provision in the AGDD geometry model for
* magnetic field information. This information was introduced in HDDS
* through a new concept called a "region" and a new tag "apply" which
* associates a volume in the geometry definition with a region. For more
* information about the meaning of regions in HDDS and the methods for
* specifying a region's magnetic field, see the documentation in the
* HDDS schema file. Geant needs to distinguish between 4 cases:
* ifield=0 : no magnetic field
* ifield=1 : general case of inhomogenous field (Runge-Kutta)
* ifield=2 : quasi-homogenous field with map (helical segments)
* ifield=3 : uniform field (helices along local z-axis)
* For detector regions with no magnetic field, the volume is created with
* ifield=0. If a field map is provided for the region, it is assigned
* ifield=1 and the field interpolator is written as a part of the output
* code. If the field is uniform then it is assigned ifield=2 unless it
* is along the z-axis, in which case it is assigned ifield=3.
*
* 3. What to do about stackX/stackY/stackZ tags
*
* In the case of Boolean tags (union/intersection/subtraction) the choice
* was easy: no support in Geant3 for these cases. For the stacks it is
* possible to construct such structures in Geant3 but it is complicated
* by the fact that the stacks do not give information about the kind or
* size of volume that should be used for the mother. Since stacks can
* be implemented easily using compositions anyway, I decided not to include
* support for them in hdds-geant.
*/
#include <xercesc/util/PlatformUtils.hpp>
#include <xercesc/util/XMLString.hpp>
#include <xercesc/util/XMLStringTokenizer.hpp>
#include <xercesc/sax/SAXParseException.hpp>
#include <xercesc/parsers/XercesDOMParser.hpp>
#include <xercesc/framework/LocalFileFormatTarget.hpp>
#include <xercesc/dom/DOM.hpp>
#include <xercesc/util/XercesDefs.hpp>
#include <xercesc/sax/ErrorHandler.hpp>
using namespace xercesc;
#include "XString.hpp"
#include "XParsers.hpp"
#include "hddsCommon.hpp"
#include <assert.h>
#include <stdlib.h>
#include <ctype.h>
#include <stdio.h>
#include <math.h>
#include <iostream>
#include <fstream>
#include <sstream>
#include <string>
#include <iomanip>
#include <vector>
#include <list>
#define APP_NAME "hddsCommon"
#define X(str) XString(str).unicode_str()
#define S(str) str.c_str()
#define NOT_USED(x) ((void)(x))
/* Refsys class:
* Stores persistent information about coordinate
* reference systems which are used to orient and place
* volumes in the simulation geometry.
*/
int Refsys::fVolumes = 0;
int Refsys::fRegions = 0;
int Refsys::fRotations = 0;
std::map<std::string,Refsys::VolIdent> Refsys::fIdentifiers;
std::vector<std::map<std::string,std::vector<int> > > Refsys::fIdentifierTable;
Refsys::Refsys() // empty constructor
: fMother(0),
fRegion(0),
fPhiOffset(0),
fRegionID(0),
fGeometryLayer(0),
fRelativeLayer(0),
fIdentifier()
{
fMOrigin[0] = fMOrigin[1] = fMOrigin[2] = 0;
fMRmatrix[0][0] = fMRmatrix[1][1] = fMRmatrix[2][2] = 1;
fMRmatrix[0][1] = fMRmatrix[1][0] = fMRmatrix[1][2] =
fMRmatrix[0][2] = fMRmatrix[2][0] = fMRmatrix[2][1] = 0;
reset();
}
Refsys::Refsys(const Refsys& src) // copy constructor
: fMother(src.fMother),
fRegion(src.fRegion),
fPhiOffset(src.fPhiOffset),
fRegionID(src.fRegionID),
fGeometryLayer(src.fGeometryLayer),
fRelativeLayer(src.fRelativeLayer),
fIdentifier(src.fIdentifier),
fPartition(src.fPartition)
{
for (int i=0; i<3; i++)
{
fMOrigin[i] = src.fMOrigin[i];
fMRmatrix[i][0] = src.fMRmatrix[i][0];
fMRmatrix[i][1] = src.fMRmatrix[i][1];
fMRmatrix[i][2] = src.fMRmatrix[i][2];
}
std::map<std::string,double>::const_iterator iter;
for (iter = src.fPar.begin(); iter != src.fPar.end(); ++iter)
{
fPar[iter->first] = iter->second;
}
reset(src);
}
Refsys& Refsys::operator=(const Refsys& src) // copy operator (deep sematics)
{
fIdentifier = src.fIdentifier;
fMother = src.fMother;
fRegion = src.fRegion;
fRegionID = src.fRegionID;
fGeometryLayer = src.fGeometryLayer;
fRelativeLayer = src.fRelativeLayer;
fPhiOffset = src.fPhiOffset;
fPartition = src.fPartition;
for (int i=0; i<3; i++)
{
fMOrigin[i] = src.fMOrigin[i];
fMRmatrix[i][0] = src.fMRmatrix[i][0];
fMRmatrix[i][1] = src.fMRmatrix[i][1];
fMRmatrix[i][2] = src.fMRmatrix[i][2];
}
std::map<std::string,double>::const_iterator iter;
for (iter = src.fPar.begin(); iter != src.fPar.end(); ++iter)
{
fPar[iter->first] = iter->second;
}
reset(src);
return *this;
}
Refsys& Refsys::reset() // reset origin, Rmatrix to null
{
fOrigin[0] = fOrigin[1] = fOrigin[2] = 0;
fRmatrix[0][0] = fRmatrix[1][1] = fRmatrix[2][2] = 1;
fRmatrix[0][1] = fRmatrix[1][0] = fRmatrix[1][2] =
fRmatrix[0][2] = fRmatrix[2][0] = fRmatrix[2][1] = 0;
fRotation = 0;
return *this;
}
Refsys& Refsys::reset(const Refsys& ref) // reset origin, Rmatrix to ref
{
for (int i = 0; i < 3; i++)
{
fOrigin[i] = ref.fOrigin[i];
fRmatrix[i][0] = ref.fRmatrix[i][0];
fRmatrix[i][1] = ref.fRmatrix[i][1];
fRmatrix[i][2] = ref.fRmatrix[i][2];
}
fRotation = ref.fRotation;
return *this;
}
Refsys& Refsys::shift(const double vector[3]) // translate origin
{
for (int i = 0; i < 3; i++)
{
fOrigin[i] += fRmatrix[i][0] * vector[0] +
fRmatrix[i][1] * vector[1] +
fRmatrix[i][2] * vector[2];
fMOrigin[i] += fMRmatrix[i][0] * vector[0] +
fMRmatrix[i][1] * vector[1] +
fMRmatrix[i][2] * vector[2];
}
return *this;
}
Refsys& Refsys::shift(const Refsys& ref) // copy origin from ref
{
fOrigin[0] = ref.fOrigin[0];
fOrigin[1] = ref.fOrigin[1];
fOrigin[2] = ref.fOrigin[2];
fMOrigin[0] = ref.fMOrigin[0];
fMOrigin[1] = ref.fMOrigin[1];
fMOrigin[2] = ref.fMOrigin[2];
return *this;
}
Refsys& Refsys::shift(const Refsys& ref,
const double vector[3]) // translate origin in ref frame
{
Refsys myRef(ref);
myRef.shift(vector);
return shift(myRef);
}
Refsys& Refsys::rotate(const double omega[3]) // rotate by vector omega (rad)
{
if ( (omega[0] != 0) || (omega[1] != 0) || (omega[2] != 0) )
{
double cosx = cos(omega[0]);
double sinx = sin(omega[0]);
double cosy = cos(omega[1]);
double siny = sin(omega[1]);
double cosz = cos(omega[2]);
double sinz = sin(omega[2]);
for (int i = 0; i < 3; i++)
{
double x[3];
double xx[3];
x[0] = fRmatrix[i][0] * cosz + fRmatrix[i][1] * sinz;
x[1] = fRmatrix[i][1] * cosz - fRmatrix[i][0] * sinz;
x[2] = fRmatrix[i][2];
xx[0] = x[0] * cosy - x[2] * siny;
xx[1] = x[1];
xx[2] = x[2] * cosy + x[0] * siny;
fRmatrix[i][0] = xx[0];
fRmatrix[i][1] = xx[1] * cosx + xx[2] * sinx;
fRmatrix[i][2] = xx[2] * cosx - xx[1] * sinx;
x[0] = fMRmatrix[i][0] * cosz + fMRmatrix[i][1] * sinz;
x[1] = fMRmatrix[i][1] * cosz - fMRmatrix[i][0] * sinz;
x[2] = fMRmatrix[i][2];
xx[0] = x[0] * cosy - x[2] * siny;
xx[1] = x[1];
xx[2] = x[2] * cosy + x[0] * siny;
fMRmatrix[i][0] = xx[0];
fMRmatrix[i][1] = xx[1] * cosx + xx[2] * sinx;
fMRmatrix[i][2] = xx[2] * cosx - xx[1] * sinx;
}
fRotation = -1;
}
return *this;
}
Refsys& Refsys::rotate(const Refsys& ref) // copy Rmatrix from ref
{
for (int i = 0; i < 3; i++)
{
fRmatrix[i][0] = ref.fRmatrix[i][0];
fRmatrix[i][1] = ref.fRmatrix[i][1];
fRmatrix[i][2] = ref.fRmatrix[i][2];
fMRmatrix[i][0] = ref.fMRmatrix[i][0];
fMRmatrix[i][1] = ref.fMRmatrix[i][1];
fMRmatrix[i][2] = ref.fMRmatrix[i][2];
}
fRotation = ref.fRotation;
return *this;
}
Refsys& Refsys::rotate(const Refsys& ref,
const double omega[3]) // rotate by omega in ref frame
{
Refsys myRef(ref);
myRef.rotate(omega);
return rotate(myRef);
}
std::vector<double>& Refsys::getRotation() const
{
std::vector<double> *angles = new std::vector<double>(3);
if (fRmatrix[2][1] == 0 && fRmatrix[2][2] == 0) {
if (fRmatrix[2][0] < 0) {
(*angles)[0] = atan2(fRmatrix[0][1],fRmatrix[1][1]);
(*angles)[1] = M_PI/2.;
(*angles)[2] = 0;
}
else {
(*angles)[0] = atan2(-fRmatrix[0][1],fRmatrix[1][1]);
(*angles)[1] = -M_PI/2.;
(*angles)[2] = 0;
}
}
else {
(*angles)[0] = atan2(fRmatrix[2][1],fRmatrix[2][2]);
(*angles)[1] = atan2(-fRmatrix[2][0],
fRmatrix[2][2]/(cos((*angles)[0])+1e-100));
(*angles)[2] = atan2(fRmatrix[1][0],fRmatrix[0][0]);
}
return *angles;
}
std::vector<double>& Refsys::getMRotation() const
{
std::vector<double> *angles = new std::vector<double>(3);
if (fMRmatrix[2][1] == 0 && fMRmatrix[2][2] == 0) {
if (fMRmatrix[2][0] < 0) {
(*angles)[0] = atan2(fMRmatrix[0][1],fMRmatrix[1][1]);
(*angles)[1] = M_PI/2.;
(*angles)[2] = 0;
}
else {
(*angles)[0] = atan2(-fMRmatrix[0][1],fMRmatrix[1][1]);
(*angles)[1] = -M_PI/2.;
(*angles)[2] = 0;
}
}
else {
(*angles)[0] = atan2(fMRmatrix[2][1],fMRmatrix[2][2]);
(*angles)[1] = atan2(-fMRmatrix[2][0],
fMRmatrix[2][2]/(cos((*angles)[0])+1e-100));
(*angles)[2] = atan2(fMRmatrix[1][0],fMRmatrix[0][0]);
}
return *angles;
}
void Refsys::clearIdentifiers()
{
fIdentifier.clear();
}
void Refsys::incrementIdentifiers()
{
std::map<std::string,Refsys::VolIdent>::iterator iter;
for (iter = fIdentifier.begin(); iter != fIdentifier.end(); ++iter)
{
iter->second.value += iter->second.step;
}
}
void Refsys::addIdentifier(XString ident, int value, int step)
{
VolIdent id;
id.value = value;
id.step = step;
fIdentifier[ident] = id;
fIdentifiers[ident] = id;
}
int Refsys::nextRotationID()
{
return ++fRotations;
}
int Refsys::nextVolumeID()
{
unsigned int ivolu = ++fVolumes;
while (fIdentifierTable.size() <= ivolu)
{
std::map<std::string,std::vector<int> > unmarked;
fIdentifierTable.push_back(unmarked);
}
return ivolu;
}
int Refsys::nextRegionID()
{
return ++fRegions;
}
/* Substance class:
* Computes and saves properties of materials that are used
* in the detector description, sometimes using directly the
* values in the hdds file and other times computing them
* based on the hdds information.
*/
Substance::Substance()
: fUniqueID(0),
fBrewList(0),
fMaterialEl(0),
fAtomicWeight(0),
fAtomicNumber(0),
fDensity(-1),
fRadLen(0),
fAbsLen(0),
fColLen(0),
fMIdEdx(0)
{}
Substance::Substance(const Substance& src)
: fUniqueID(src.fUniqueID),
fMaterialEl(src.fMaterialEl),
fAtomicWeight(src.fAtomicWeight),
fAtomicNumber(src.fAtomicNumber),
fDensity(src.fDensity),
fRadLen(src.fRadLen),
fAbsLen(src.fAbsLen),
fColLen(src.fColLen),
fMIdEdx(src.fMIdEdx)
{
std::list<Brew>::const_iterator iter;
for (iter = src.fBrewList.begin();
iter != src.fBrewList.end();
++iter)
{
Brew formula = *iter;
formula.sub = new Substance(*iter->sub);
fBrewList.push_back(formula);
}
}
Substance::Substance(DOMElement* elem)
: fUniqueID(0),
fBrewList(0),
fMaterialEl(elem),
fAtomicWeight(0),
fAtomicNumber(0),
fDensity(-1),
fRadLen(0),
fAbsLen(0),
fColLen(0),
fMIdEdx(0)
{
XString aS(fMaterialEl->getAttribute(X("a")));
fAtomicWeight = atof(S(aS));
XString zS(fMaterialEl->getAttribute(X("z")));
fAtomicNumber = atof(S(zS));
double wfactSum = 0;
DOMNode* cont;
for (cont = fMaterialEl->getFirstChild(); cont != 0;
cont = cont->getNextSibling())
{
if (cont->getNodeType() == DOMNode::ELEMENT_NODE)
{
DOMElement* contEl = (DOMElement*)cont;
XString tagS(contEl->getTagName());
if (tagS == "real")
{
Units unit;
unit.getConversions(contEl);
XString nameS(contEl->getAttribute(X("name")));
XString valueS(contEl->getAttribute(X("value")));
if (nameS == "density")
{
fDensity = atof(S(valueS)) /(unit.g/unit.cm3);
}
else if (nameS == "radlen")
{
fRadLen = atof(S(valueS)) /unit.cm;
}
else if (nameS == "abslen")
{
fAbsLen = atof(S(valueS)) /unit.cm;
}
else if (nameS == "collen")
{
fColLen = atof(S(valueS)) /unit.cm;
}
else if (nameS == "dedx")
{
fMIdEdx = atof(S(valueS)) /(unit.MeV/unit.cm);
}
}
else if (tagS == "addmaterial")
{
XString matS(contEl->getAttribute(X("material")));
DOMDocument* document = fMaterialEl->getOwnerDocument();
DOMElement* targEl = document->getElementById(X(matS));
XString typeS(targEl->getTagName());
Substance::Brew formula;
formula.sub = new Substance(targEl);
formula.natoms = 0;
formula.wfact = 0;
DOMNode* mix;
for (mix = contEl->getFirstChild(); mix != 0;
mix = mix->getNextSibling())
{
if (mix->getNodeType() == DOMNode::ELEMENT_NODE)
{
DOMElement* mixEl = (DOMElement*)mix;
XString mixS(mixEl->getTagName());
if (mixS == "natoms")
{
if (typeS != "element")
{
std::cerr << APP_NAME
<< " - error processing composite "
<< S(matS) << std::endl
<< "natoms can only be specified for elements."
<< std::endl;
exit(1);
}
XString nS(mixEl->getAttribute(X("n")));
formula.natoms = atoi(S(nS));
formula.wfact = formula.natoms * formula.sub->getAtomicWeight();
}
else if (mixS == "fractionmass")
{
XString fS(mixEl->getAttribute(X("fraction")));
formula.wfact = atof(S(fS));
}
}
}
fBrewList.push_back(formula);
wfactSum += formula.wfact;
}
}
}
double aSum=0; double aNorm=0;
double zSum=0; double zNorm=0;
double densitySum=0; double densityNorm=0;
double radlenSum=0; double radlenNorm=0;
double abslenSum=0; double abslenNorm=0;
double collenSum=0; double collenNorm=0;
double dedxSum=0; double dedxNorm=0;
std::list<Substance::Brew>::iterator iter;
for (iter = fBrewList.begin(); iter != fBrewList.end(); ++iter)
{
iter->wfact /= wfactSum;
double weight = (iter->natoms > 0)?
iter->natoms * iter->sub->fAtomicWeight : iter->wfact;
aSum += weight*iter->sub->fAtomicWeight;
aNorm += weight;
zSum += weight*iter->sub->fAtomicNumber;
zNorm += weight;
densitySum += weight/iter->sub->fDensity;
densityNorm += weight;
weight /= iter->sub->fDensity;
radlenSum += weight/(iter->sub->fRadLen + 1e-99);
radlenNorm += weight;
abslenSum += weight/(iter->sub->fAbsLen + 1e-99);
abslenNorm += weight;
collenSum += weight/(iter->sub->fColLen + 1e-99);
collenNorm += weight;
dedxSum += weight*iter->sub->fMIdEdx;
dedxNorm += weight;
}
if (fAtomicWeight == 0 && aNorm > 0)
{
fAtomicWeight = aSum/aNorm;
}
if (fAtomicNumber == 0 && zNorm > 0)
{
fAtomicNumber = zSum/zNorm;
}
if (fDensity <= 0 && densitySum > 0)
{
fDensity = densityNorm/densitySum;
}
if (fRadLen == 0 && radlenSum > 0)
{
fRadLen = radlenNorm/radlenSum;
}
if (fAbsLen == 0 && abslenSum > 0)
{
fAbsLen = abslenNorm/abslenSum;
}
if (fColLen == 0 && collenSum > 0)
{
fColLen = collenNorm/collenSum;
}
if (fMIdEdx == 0 && dedxNorm > 0)
{
fMIdEdx = dedxSum/dedxNorm;
}
if (fDensity < 0)
{
XString tagS(fMaterialEl->getTagName());
XString nameS(fMaterialEl->getAttribute(X("name")));
std::cerr
<< APP_NAME << " error: " << S(tagS) << " " << S(nameS)
<< ", atomic number " << fAtomicNumber
<< ", atomic weight " << fAtomicWeight
<< " is missing a density specification." << std::endl;
exit(1);
}
}
Substance::~Substance()
{
std::list<Brew>::iterator iter;
for (iter = fBrewList.begin(); iter != fBrewList.end(); ++iter)
{
delete iter->sub;
}
}
Substance& Substance::operator=(const Substance& src)
{
fMaterialEl = src.fMaterialEl;
fAtomicWeight = src.fAtomicWeight;
fAtomicNumber = src.fAtomicNumber;
fDensity = src.fDensity;
fRadLen = src.fRadLen;
fAbsLen = src.fAbsLen;
fColLen = src.fColLen;
fMIdEdx = src.fMIdEdx;
fUniqueID = src.fUniqueID;
fBrewList = src.fBrewList;
std::list<Brew>::iterator iter;
for (iter = fBrewList.begin();
iter != fBrewList.end();
++iter)
{
iter->sub = new Substance(*iter->sub);
}
return *this;
}
double Substance::getAtomicWeight()
{
return fAtomicWeight;
}
double Substance::getAtomicNumber()
{
return fAtomicNumber;
}
double Substance::getDensity()
{
return fDensity;
}
double Substance::getRadLength()
{
return fRadLen;
}
double Substance::getAbsLength()
{
return fAbsLen;
}
double Substance::getColLength()
{
return fColLen;
}
double Substance::getMIdEdx()
{
return fMIdEdx;
}
XString Substance::getName()
{
return XString(fMaterialEl->getAttribute(X("name")));
}
XString Substance::getSymbol()
{
return XString(fMaterialEl->getAttribute(X("symbol")));
}
DOMElement* Substance::getDOMElement()
{
return fMaterialEl;
}
/* Units class:
* Provides conversion constants for convenient extraction of
* dimensioned parameters in the hdds file.
*/
Units::Units()
{
set_1s(1.);
set_1cm(1.);
set_1rad(1.);
set_1MeV(1.);
set_1g(1.);
set_1cm2(1.);
set_1cm3(1.);
set_1G(1.);
percent = 0.01;
}
Units::Units(const Units& u)
{
set_1s(u.s);
set_1cm(u.cm);
set_1rad(u.rad);
set_1MeV(u.MeV);
set_1g(u.g);
set_1cm2(u.cm2);
set_1cm3(u.cm3);
set_1G(u.G);
percent = 0.01;
}
void Units::set_1s(double tu)
{
s=tu;
ns=s*1e-9; ms=s*1e-3;
min=s/60; hr=min/60; days=hr/24; weeks=days/7;
}
void Units::set_1cm(double lu)
{
cm=lu;
m=cm*1e2; mm=m*1e-3; um=m*1e-6; nm=m*1e-9; km=m*1e3;
in=cm*2.54; ft=in*12; miles=ft*5280; mils=in*1e-3;
}
void Units::set_1rad(double au)
{
rad=au;
mrad=rad*1e-3; urad=rad*1e-6;
deg=rad*M_PI/180; arcmin=deg/60; arcsec=arcmin/60;
}
void Units::set_1deg(double au)
{
deg=au;
arcmin=deg/60; arcsec=arcmin/60;
rad=deg*180/M_PI; mrad=rad*1e-3; urad=rad*1e-6;
}
void Units::set_1MeV(double eu)
{
MeV=eu;
eV=MeV*1e-6; KeV=eV*1e3; GeV=eV*1e9; TeV=eV*1e12;
}
void Units::set_1g(double mu)
{
g=mu;
kg=g*1e3; mg=g*1e-3;
}
void Units::set_1cm2(double l2u)
{
cm2=l2u;
m2=cm2*1e4; mm2=cm2*1e-2;
b=cm2*1e-24; mb=b*1e-3; ub=b*1e-6; nb=b*1e-9; pb=b*1e-12;
}
void Units::set_1cm3(double l3u)
{
cm3=l3u;
ml=cm3; l=ml*1e3;
}
void Units::set_1G(double bfu)
{
G=bfu;
kG=G*1e3; Tesla=kG*10;
}
void Units::getConversions(DOMElement* el)
{
XString unitlS(el->getAttribute(X("unit_length")));
if (unitlS.size() == 0)
{
;
}
else if (unitlS == "mm")
{
set_1cm(10);
}
else if (unitlS == "cm")
{
set_1cm(1);
}
else if (unitlS == "m")
{
set_1cm(0.01);
}
else if (unitlS == "km")
{
set_1cm(1e-5);
}
else if (unitlS == "um")
{
set_1cm(1e6);
}
else if (unitlS == "nm")
{
set_1cm(1e9);
}
else if (unitlS == "in")
{
set_1cm(1/2.54);
}
else if (unitlS == "ft")
{
set_1cm(1/(12*2.54));
}
else if (unitlS == "mil")
{
set_1cm(1000/2.54);
}
else
{
XString tagS(el->getTagName());
std::cerr
<< APP_NAME << " error: unknown length unit " << S(unitlS)
<< " on tag " << S(tagS) << std::endl;
exit(1);
}
XString unitaS = el->getAttribute(X("unit_angle"));
if (unitaS.size() == 0)
{
;
}
else if (unitaS == "deg")
{
set_1deg(1);
}
else if (unitaS == "rad")
{
set_1rad(1);
}
else if (unitaS == "mrad")
{
set_1rad(1e3);
}
else
{
XString tagS(el->getTagName());
std::cerr
<< APP_NAME << " error: unknown angle unit " << S(unitaS)
<< " on volume " << S(tagS) << std::endl;
exit(1);
}
XString unitS = el->getAttribute(X("unit"));
if (unitS.size() == 0)
{
;
}
else if (unitS == "mm")
{
set_1cm(10);
}
else if (unitS == "cm")
{
set_1cm(1);
}
else if (unitS == "m")
{
set_1cm(0.01);
}
else if (unitlS == "km")
{
set_1cm(1e-5);
}
else if (unitlS == "um")
{
set_1cm(1e6);
}
else if (unitlS == "nm")
{
set_1cm(1e9);
}
else if (unitlS == "in")
{
set_1cm(1/2.54);
}
else if (unitlS == "ft")
{
set_1cm(1/(12*2.54));
}
else if (unitlS == "mil")
{
set_1cm(1000/2.54);
}
else if (unitaS == "deg")
{
set_1deg(1);
}
else if (unitaS == "rad")
{
set_1rad(1);
}
else if (unitaS == "mrad")
{
set_1rad(1e3);
}
else if (unitS == "eV")
{
set_1MeV(1e6);
}
else if (unitS == "KeV")
{
set_1MeV(1e3);
}
else if (unitS == "MeV")
{
set_1MeV(1);
}
else if (unitS == "GeV")
{
set_1MeV(1e-3);
}
else if (unitS == "g/cm^2")
{
set_1g(1);
set_1cm2(1);
}
else if (unitS == "g/cm^3")
{
set_1g(1);
set_1cm3(1);
}
else if (unitS == "MeV/g/cm^2")
{
set_1MeV(1);
set_1g(1);
set_1cm2(1);
}
else if (unitS == "Tesla" || unitS == "T")
{
set_1G(1e-4);
}
else if (unitS == "kG" || unitS == "kGs")
{
set_1G(1e-3);
}
else if (unitS == "G" || unitS == "Gs")
{
set_1G(1);
}
else if (unitS == "percent")
{
;
}
else if (unitS == "none")
{
;