StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
StiVMCToolKit.cxx
1 //#define __NOVMC__
2 #include "StiVMCToolKit.h"
3 #ifdef __ROOT__
4 #include "StMaker.h"
5 #endif
6 #if !defined(__CINT__) || defined(__MAKECINT__)
7 #include <string.h>
8 #include <assert.h>
9 #include "Riostream.h"
10 #include <stdio.h>
11 #include "TSystem.h"
12 #include "TMath.h"
13 #include "TFile.h"
14 #include "StMessMgr.h"
15 #else
16 #define BIT(n) (1 << (n))
17 #endif
18 static Int_t __addelement__(Int_t NElem,const TGeoMaterial *mat, Elem_t *ElementList) {
19  return StiVMCToolKit::Add2ElementList(NElem, mat, ElementList);
20 }
21 
22 struct Elem_t {
23  Int_t index;
24  Double_t W;
25  Double_t A;
26  Double_t Z;
27  Elem_t() : index(0),W(0),A(0),Z(0) {}
28 };
29 struct ElemV_t {
30  const Char_t *name;
31  Double_t A;
32  Double_t Z;
33  Double_t V; // by volume
34  Double_t W;
35 };
36 static Int_t m_Debug = 0;
37 static const Int_t NoElemMax = 10000;
38 static VolumeMap_t VolumesToBeAveraged[] = {
39  {"PIPE","the STAR beam pipe mother volume","HALL_1/CAVE_1/PIPE_1-2/*","",""},
40  {"PIPC","the Central Beam PIPe Volum","HALL_1/CAVE_1/PIPE_1-2/PIPC_1/*","",""},
41  {"PVAC","the Vacuum Volume of Be section of pipe","HALL_1/CAVE_1/PIPE_1-2/PIPC_1/PVAC_1","",""},
42  {"PIPO","Steel pipe from Be to 1st flanges","HALL_1/CAVE_1/PIPE_1-2/PIPO_1/PVAO_1","",""},
43  {"PVAO","its cavity","HALL_1/CAVE_1/PIPE_1-2/PIPO_1/PVAO_1","",""},
44  // Pixel
45  {"PXBX","Extrenal Berillium tube","HALL_1/CAVE_1/PXBX_1","",""},
46  // SVT
47  {"SOUM", "Outer shileding structure","HALL_1/CAVE_1/SVTT_1/SOUM_1/*","",""},
48  {"SXRL", "Circular water feeds","HALL_1/CAVE_1/SVTT_1/SXRL_1-2/*","",""},
49  {"SXR1", "Circular water feeds","HALL_1/CAVE_1/SVTT_1/SXR1_3-4/*","",""},
50  {"SXR2", "Circular water feeds","HALL_1/CAVE_1/SVTT_1/SXR2_5-6/*","",""},
51  {"SCBM", "Mother of All Cables","HALL_1/CAVE_1/SVTT_1/SCBM_1-2/*","",""},
52  {"SCBL", "The bundles of cables connecting PCBs with the transition boards","HALL_1/CAVE_1/SVTT_1/SCBM_1-2/SCBL_1","",""},
53  {"SCB1", "The bundles of cables connecting PCBs with the transition boards","HALL_1/CAVE_1/SVTT_1/SCBM_1-2/SCB1_2","",""},
54  {"SCB2", "The bundles of cables connecting PCBs with the transition boards","HALL_1/CAVE_1/SVTT_1/SCBM_1-2/SCB2_3","",""},
55  {"SCB3", "The bundles of cables connecting PCBs with the transition boards","HALL_1/CAVE_1/SVTT_1/SCBM_1-2/SCB3_4","",""},
56  {"SFED", "bundles of water pipes","HALL_1/CAVE_1/SVTT_1/SCBM_1-2/SFED_1","",""},
57  {"SFE1", "bundles of water pipes","HALL_1/CAVE_1/SVTT_1/SCBM_1-2/SFE1_2","",""},
58  {"SFE2", "bundles of water pipes","HALL_1/CAVE_1/SVTT_1/SCBM_1-2/SFE2_3","",""},
59  {"SPLS", "plastic of the water pipes","HALL_1/CAVE_1/SVTT_1/SCBM_1-2/SPLS_1","",""},
60  {"SPL1", "plastic of the water pipes","HALL_1/CAVE_1/SVTT_1/SCBM_1-2/SPL1_2","",""},
61  {"SPL2", "plastic of the water pipes","HALL_1/CAVE_1/SVTT_1/SCBM_1-2/SPL2_3","",""},
62  {"SALM", "aluminum shield mesh","HALL_1/CAVE_1/SVTT_1/SALM_1-2","",""},
63  {"SOSH", "SVT outer shield","HALL_1/CAVE_1/SVTT_1/SOSH_1","",""},
64  {"SISH", "SVT inner shield","HALL_1/CAVE_1/SVTT_1/SISH_1","",""},
65  {"SLYD", "layer mother","HALL_1/CAVE_1/SVTT_1/SLYD_1/*","",""},
66  {"SLY1", "layer mother","HALL_1/CAVE_1/SVTT_1/SLY1_2/*","",""},
67  {"SLY2", "layer mother","HALL_1/CAVE_1/SVTT_1/SLY2_3/*","",""},
68  {"SLY3", "layer mother","HALL_1/CAVE_1/SVTT_1/SLY3_4/*","",""},
69  {"SLY4", "layer mother","HALL_1/CAVE_1/SVTT_1/SLY4_5/*","",""},
70  {"SLY5", "layer mother","HALL_1/CAVE_1/SVTT_1/SLY5_6/*","",""},
71  {"SVTD", "an active wafer volume","HALL_1/CAVE_1/SVTT_1/SLY*/SLS*/SLD*/STL*/STS*/SVTD_1","svt","SVTD"}, // <+++
72  {"SROD", "Support rod","HALL_1/CAVE_1/SVTT_1/SROD_1-2","",""},
73  {"SBSP", "Beampipe support mother","HALL_1/CAVE_1/SVTT_1/SBSP_1-2","",""},
74  // {"SCON", "Support cone mother","HALL_1/CAVE_1/SVTT_1/SCON_1-2/*","",""},
75  {"SBWC", "water manifold to support cone bracket mother","HALL_1/CAVE_1/SVTT_1/SBWC_1-2/*","",""},
76  {"SWMM", "water manifold mother","HALL_1/CAVE_1/SVTT_1/SWMM_1-2/*","",""},
77  {"SIES", "Volume to hold inner endring screws","HALL_1/CAVE_1/SVTT_1/SIES_1-2/*","",""},
78  {"SOES", "Volume to hold outer endring screws","HALL_1/CAVE_1/SVTT_1/SOES_1-2/*","",""},
79  {"SBRG", "Bracket joining the end rungs","HALL_1/CAVE_1/SVTT_1/SBRG_1-2/*","",""},
80  {"SOER", "outer end ring","HALL_1/CAVE_1/SVTT_1/SOER_1-2/*","",""},
81  {"SIRT", "inner end ring tube piece ","HALL_1/CAVE_1/SVTT_1/SIRT_1-2","",""},
82  {"SIRP", "inner end ring polygon piece ","HALL_1/CAVE_1/SVTT_1/SIRP_1-2","",""},
83  {"STAC", "twinax cable approximation, copper","HALL_1/CAVE_1/SVTT_1/SCON_1/STAC_1-2","",""},
84  // SSD
85  // {"SFMO", "the mother of all Silicon Strip Detector volumes","HALL_1/CAVE_1/SVTT_1/SFMO_1","",""},
86  {"SCMP","SSD mounting plate inserted in the cone","HALL_1/CAVE_1/SVTT_1/SFMO_1/SCMP_1-8","",""},
87  {"SCVM","SSD V-shape mouting piece","HALL_1/CAVE_1/SVTT_1/SFMO_1/SCVM_1-8/*","",""},
88  {"SSLT","the linking (sector to the cone) tube","HALL_1/CAVE_1/SVTT_1/SFMO_1/SSLT_1-8","",""},
89  {"SSLB","the linking (sector to the cone)","HALL_1/CAVE_1/SVTT_1/SFMO_1/SSLB_1-8","",""},
90  {"SSRS","the side of the small rib","HALL_1/CAVE_1/SVTT_1/SFMO_1/SSRS_1-4","",""},
91  {"SSRT","the top of the side rib","HALL_1/CAVE_1/SVTT_1/SFMO_1/SSRT_1-4","",""},
92  {"SSSS","Side parts of the small sectors","HALL_1/CAVE_1/SVTT_1/SFMO_1/SSSS_1-4","",""},
93  {"SSST","Top parts of the small sectors","HALL_1/CAVE_1/SVTT_1/SFMO_1/SSST_1-4","",""},
94  // {"SFLM","the mother of the ladder","HALL_1/CAVE_1/SVTT_1/SFMO_1/SFLM_1-20/*","",""},
95  {"SFSM","the structure mother volume","HALL_1/CAVE_1/SVTT_1/SFMO_1/SFLM_1-20/SFSM_1/*","",""},
96  {"SFDM","the detectors and adcs mother volume","HALL_1/CAVE_1/SVTT_1/SFMO_1/SFLM_1-20/SFDM_1/*","",""},
97  {"SFSD","the strip detector", "HALL_1/CAVE_1/SVTT_1/SFMO_1/SFLM_1-20/SFDM_1/SFSW_1-16/SFSD_1","ssd",""},// <+++
98  // TPC
99  // {"TPCE","the TPC system in STAR","HALL_1/CAVE_1/TPCE_1","",""},
100  {"TPCW","the TPC supporting endcap Wheel","HALL_1/CAVE_1/TPCE_1/TPCW_1-2/*","",""},
101  {"TPEA","one endcap placed in TPC","HALL_1/CAVE_1/TPCE_1/TPEA_1-2/*","",""},
102  {"TPCM","the Central Membrane placed in TPC","HALL_1/CAVE_1/TPCE_1/TPCM_1","",""},
103  {"TOFC","defines outer field cage - fill it with insulating gas already","HALL_1/CAVE_1/TPCE_1/TOFC_1/*","",""},
104  {"TIFC","defines the Inner Field Cage placed in TPC","HALL_1/CAVE_1/TPCE_1/TIFC_1/*","",""},
105  {"TPGV","the Gas Volume placed in TPC","HALL_1/CAVE_1/TPCE_1/TPGV_1-2/*","",""},
106  {"TPSS","a division of gas volume corresponding to a supersectors","HALL_1/CAVE_1/TPCE_1/TPGV_1-2/TPSS_1-12/*","",""},
107  {"TPAD","(inner) real padrow with dimensions defined at positioning time","HALL_1/CAVE_1/TPCE_1/TPGV_1-2/TPSS_1-12/TPAD_1-39","tpc",""},// <+++
108  {"TPA1","(outer) real padrow with dimensions defined at positioning time","HALL_1/CAVE_1/TPCE_1/TPGV_1-2/TPSS_1-12/TPA1_40-73","tpc",""}
109 };
110 static Int_t NofVolToBEAveraged = sizeof(VolumesToBeAveraged)/sizeof(VolumeMap_t);
111 static VolumeMap_t TopVolumes[] = {// Mother volume and sensitive detectors
112  {"BBCM","Beam Beam Counter Modules Geometry","","",""},
113  {"BTOF","the whole CTF system envelope","","",""},
114  {"CALB","the geometry of the Barrel EM Calorimeter","","",""},
115  {"ECAL","the EM EndCap Calorimeter GEOmetry","","",""},
116  {"FBOX","one Pb-Glass fpd detector","","",""},
117  {"FBO1","an other Pb-Glass fpd detector","","",""},
118  {"FTMO","the mother of the single FTPC RO barrel","","",""},
119  {"IBEM","the IBeam structure beneath the Bell reducer cone","","",""},
120  {"MAGP","the geometry of the STAR magnet","","",""},
121  {"PHMD","the Photon Multiplicity Detector","","",""},
122  {"SUPO","the geometry of the Forward TPC supports in STAR","","",""},
123 
124  {"SVTT","the SVT geometry for STAR","","",""},
125 
126  {"SFMO","the mother of all Silicon Strip Detector volumes (inside of SVTT)","","",""},
127 
128  {"FTPC","the geometry of the Forward TPC in STAR (inside of SVTT)","","",""},
129 
130  {"UPST","the geometry of the UPSTREAM AreA.","","",""},
131  {"VPDD","the Pseudo Vertex Position Detector of STAR","","",""},
132  {"ZCAL","the geometry of the Zero deg. Quartz Calorimeter","","",""}
133 };
134 static Int_t nTopVol = sizeof(TopVolumes)/sizeof(VolumeMap_t);
135  /*
136  "Air" 4 1.214e-3 #. Huhtinen, Air 18 degr.C and 58% humidity
137  "Nitrogen" 74.94 #. Weight fraction
138  "Oxygen" 23.69 #. Weight fraction
139  "Argon" 1.29 #. Weight fraction
140  "Hydrogen" 0.08 #. Weight fraction
141  */
142 static ElemV_t Air[] = {
143  {"N_2", 14.00674, 7, 2*78.084, 0},
144  {"O_2", 15.9994 , 8, 2*(20.946 + 0.033), 0},// {"CO_2", 0, 0, 2*0.033, 0}, 2*(20.946 + 0.033), 0},
145  {"C", 12.011, 6, 0.033, 0},
146  {"Ar", 39.948, 18, 0.934, 0}
147 };
148 static Int_t NAir = 0;
149 static ElemV_t P10[] = {
150  {"Ar", 40,18,0.95744680,0},
151  {"C" , 12, 6,0.03191489,0},
152  {"H" , 1, 1,0.01063830,0}
153 };
154 static Int_t NP10 = 3;
155 // {"Ne", 20.1797, 10, 18.18e-6, 0},
156 // {"He", 4.002602, 2, 5.24e-6, 0},
157 // {"Kr", 0., 0, 1.14e-6, 0}
158 // {"Xe", 0., 0, 0.087e-6, 0}
159 // {H_2", 0., 0, 0.5e-6}
160 // {"CH_4",0., 0, 2e-6}
161 // {N_2 O", 0., 0, 0.5-e6}
162 //________________________________________________________________________________
163 void StiVMCToolKit::SetDebug(Int_t m) {m_Debug = m;}
164 //________________________________________________________________________________
165 Int_t StiVMCToolKit::Debug() {return m_Debug;}
166 //________________________________________________________________________________
167 void StiVMCToolKit::PrintShape(TGeoShape *shape) {
168  TGeoBBox *box = 0, *Box = 0;
169  TGeoTrd1 *trd1 = 0;
170  TGeoTrd2 *trd2 = 0;
171  TGeoTube *tube = 0;
172  TGeoTubeSeg *tubs = 0;
173  TGeoPcon *pcon = 0;
174  TGeoPgon *pgon = 0;
175  TGeoCone *cone = 0;
176  TGeoConeSeg *cons = 0;
177  TGeoArb8 *arb8 = 0;
178  TGeoEltu *eltu = 0;
179  Double_t *XY;
180  Double_t dZ;
181  // Double_t paramsBB[3];
182  Double_t paramsBC[4];
183  Int_t i, j;
184  Int_t Nz = 0;
185  shape->GetBoundingCylinder(paramsBC);
186  shape->ComputeBBox();
187  Box = (TGeoBBox *) shape;
188  for (Int_t bit = 24; bit >= 9; bit--) {//cout << bit << "\t";
189  if (shape->TestShapeBit(BIT(bit))) {
190  switch (BIT(bit)) {
191  case TGeoShape::kGeoBox:
192  box = (TGeoBBox *) shape;
193  cout << "Box \tdX\t" << box->GetDX() << "\tdY\t" << box->GetDY() << "\tdZ\t" << box->GetDZ()
194  << endl;
195  break;
196  case TGeoShape::kGeoTrd1:
197  trd1 = (TGeoTrd1 *) shape;
198  cout << "Trd1\tdX1\t" << trd1->GetDx1() << "\tdX2\t" << trd1->GetDx2()
199  << "\tdY\t" << trd1->GetDy() << "\tdZ\t" << trd1->GetDz()
200  << endl;
201  break;
202  case TGeoShape::kGeoTrd2:
203  trd2 = (TGeoTrd2 *) shape;
204  cout << "Trd2\tdX1\t" << trd2->GetDx1() << "\tdX2\t" << trd2->GetDx2()
205  << "\tdY1\t" << trd2->GetDy1() << "\tdY2\t" << trd2->GetDy2()
206  << "\tdZ\t" << trd2->GetDz() << endl;
207  break;
208  case TGeoShape::kGeoTubeSeg:
209  tubs = (TGeoTubeSeg *) shape;
210  cout << "Tubs\tRmin\t" << tubs->GetRmin() << "\tRmax\t" << tubs->GetRmax() << "\tdZ\t" << tubs->GetDz()
211  << "\tPhi1\t" << tubs->GetPhi1() << "\tPhi2\t" << tubs->GetPhi2()
212  << endl;
213  break;
214  case TGeoShape::kGeoTube:
215  tube = (TGeoTube *) shape;
216  cout << "Tube\tRmin\t" << tube->GetRmin() << "\tRmax\t" << tube->GetRmax() << "\tdZ\t" << tube->GetDz()
217  << endl;
218  break;
219  case TGeoShape::kGeoPcon:
220  pcon = (TGeoPcon *) shape;
221  Nz = pcon->GetNz();
222  cout << "Pcon"
223  << "\tPhi1\t" << pcon->GetPhi1() << "\tDphi\t" << pcon->GetDphi() << "\tNz\t" << Nz << endl;
224  for (i = 0; i < Nz; i++) {
225  cout << i << "\tZ\t" << pcon->GetZ(i) << "\tRmin\t" << pcon->GetRmin(i) << "\tRmax\t" << pcon->GetRmax(i) << endl;
226  }
227  cout << endl;
228  break;
229  case TGeoShape::kGeoPgon:
230  pgon = (TGeoPgon *) shape;
231  Nz = pgon->GetNz();
232  // pcon = (TGeoPcon *) shape;
233  cout << "Pgon\tNedges\t" << pgon->GetNedges()
234  << "\tPhi1\t" << pgon->GetPhi1() << "\tDphi\t" << pgon->GetDphi() << "\tNz\t" <<Nz << endl;
235  for (i = 0; i <Nz; i++) {
236  cout << i << "\tZ\t" << pgon->GetZ(i) << "\tRmin\t" << pgon->GetRmin(i) << "\tRmax\t" << pgon->GetRmax(i) << endl;
237  }
238  cout << endl;
239  break;
240  case TGeoShape::kGeoCone:
241  cone = (TGeoCone *) shape;
242  cout << "Cone\tdZ\t" << cone->GetDz()
243  << "\tRmin1\t" << cone->GetRmin1() << "\tRmax1\t" << cone->GetRmax1()
244  << "\tRmin2\t" << cone->GetRmin2() << "\tRmax2\t" << cone->GetRmax2()
245  << endl;
246  break;
247  case TGeoShape::kGeoConeSeg:
248  cons = (TGeoConeSeg *) shape;
249  cout << "Cons\tdZ\t" << cons->GetDz()
250  << "\tPhi1\t" << cons->GetPhi1() << "\tPhi2\t" << cons->GetPhi2()
251  << "\tRmin1\t" << cons->GetRmin1() << "\tRmax1\t" << cons->GetRmax1()
252  << "\tRmin2\t" << cons->GetRmin2() << "\tRmax2\t" << cons->GetRmax2()
253  << endl;
254  break;
255  case TGeoShape::kGeoArb8:
256  case TGeoShape::kGeoTrap:
257  arb8 = (TGeoArb8 *) shape;
258  XY = arb8->GetVertices();
259  dZ = arb8->GetDz();
260  for (j = 0; j < 2; j++) {
261  cout << "Trap/Arb8\tdZ\t";
262  if (j == 0) cout << -dZ;
263  else cout << dZ;
264  for (i = 4*j; i < 4*j+4; i++)
265  cout << "\t(" << XY[2*i] << "," << XY[2*i+1] << ")";
266  cout << endl;
267  }
268  break;
269  case TGeoShape::kGeoEltu:
270  eltu = (TGeoEltu *) shape;
271  cout << "Eltu\tdZ\t" << eltu->GetDz()
272  << "\tA\t" << eltu->GetA() << "\tB\t" << eltu->GetB()
273  << endl;
274  break;
275  case TGeoShape::kGeoTorus:
276  case TGeoShape::kGeoPara:
277  case TGeoShape::kGeoSph:
278  case TGeoShape::kGeoCtub:
279  default:
280  cout << bit << "\t has not yet implemented for " << shape->GetName() << endl;
281  break;
282  }
283  break;
284  }
285  }
286 }
287 //________________________________________________________________________________
288 Double_t StiVMCToolKit::GetShapeVolume(TGeoShape *shape) {
289  TGeoBBox *box = 0, *Box = 0;
290  TGeoTrd1 *trd1 = 0;
291  TGeoTrd2 *trd2 = 0;
292  TGeoTube *tube = 0;
293  TGeoTubeSeg *tubs = 0;
294  TGeoPcon *pcon = 0;
295  TGeoPgon *pgon = 0;
296  TGeoCone *cone = 0;
297  TGeoConeSeg *cons = 0;
298  TGeoEltu *eltu = 0;
299  Double_t volume = 0;
300  Double_t paramsBC[4];
301  Double_t volBB = 0;
302  Double_t volBC = 0;
303  Int_t i;
304  Int_t Nz = 0;
305  shape->GetBoundingCylinder(paramsBC);
306  shape->ComputeBBox();
307  Box = (TGeoBBox *) shape;
308  volBB = 8*Box->GetDX()*Box->GetDY()*Box->GetDZ();
309  volBC = 2*TMath::Pi()*(paramsBC[1] - paramsBC[0])*Box->GetDZ();;
310  for (Int_t bit = 24; bit >= 9; bit--) {//cout << bit << "\t";
311  if (shape->TestShapeBit(BIT(bit))) {
312  Int_t kbit = BIT(bit);
313  switch (kbit) {
314  case TGeoShape::kGeoBox:
315  box = (TGeoBBox *) shape;
316  volume = 8*box->GetDX()*box->GetDY()*box->GetDZ();
317  break;
318  case TGeoShape::kGeoTrd1:
319  trd1 = (TGeoTrd1 *) shape;
320  volume = 4*(trd1->GetDx1()+trd1->GetDx2())*trd1->GetDy()*trd1->GetDz();
321  break;
322  case TGeoShape::kGeoTrd2:
323  trd2 = (TGeoTrd2 *) shape;
324  volume = 2*(trd2->GetDx1() + trd2->GetDx2())*(trd2->GetDy1() + trd2->GetDy2())*trd2->GetDz();
325  break;
326  case TGeoShape::kGeoTubeSeg:
327  case TGeoShape::kGeoTube:
328  tube = (TGeoTube *) shape;
329  volume = 2*TMath::Pi()*(tube->GetRmax()* tube->GetRmax() - tube->GetRmin()*tube->GetRmin())*
330  tube->GetDz();
331  if (kbit == TGeoShape::kGeoTubeSeg) {
332  tubs = (TGeoTubeSeg *) shape;
333  volume *= TMath::Abs(tubs->GetPhi2()-tubs->GetPhi1())/360.;
334  }
335  break;
336  case TGeoShape::kGeoPcon:
337  case TGeoShape::kGeoPgon:
338  pcon = (TGeoPcon *) shape;
339  Nz = pcon->GetNz();
340  volume = 0;
341  for (i = 1; i < Nz; i++) {
342  volume +=
343  ((pcon->GetRmax(i)*pcon->GetRmax(i) + pcon->GetRmax(i-1)*pcon->GetRmax(i-1) + pcon->GetRmax(i)*pcon->GetRmax(i-1)) -
344  (pcon->GetRmin(i)*pcon->GetRmin(i) + pcon->GetRmin(i-1)*pcon->GetRmin(i-1) + pcon->GetRmin(i)*pcon->GetRmin(i-1)))*
345  TMath::Abs((pcon->GetZ(i) - pcon->GetZ(i-1)));
346  }
347  if (kbit == TGeoShape::kGeoPgon) {
348  pgon = (TGeoPgon *) shape;
349  volume *= TMath::Tan(TMath::Pi()/180.*TMath::Abs(pgon->GetDphi()/pgon->GetNedges()/2.))/3.*pgon->GetNedges();
350  }
351  else {
352  volume *= TMath::Pi()/3.*TMath::Abs(pcon->GetDphi())/360.;
353  }
354  break;
355  case TGeoShape::kGeoCone:
356  case TGeoShape::kGeoConeSeg:
357  cone = (TGeoCone *) shape;
358  volume = TMath::Pi()*2./3.*
359  ((cone->GetRmax2()*cone->GetRmax2() + cone->GetRmax1()*cone->GetRmax1() + cone->GetRmax2()*cone->GetRmax1()) -
360  (cone->GetRmin2()*cone->GetRmin2() + cone->GetRmin1()*cone->GetRmin1() + cone->GetRmin2()*cone->GetRmin1()))*
361  cone->GetDz();
362  if (kbit == TGeoShape::kGeoConeSeg) {
363  cons = (TGeoConeSeg *) shape;
364  volume *= TMath::Abs(cons->GetPhi2() - cons->GetPhi1())/360.;
365  }
366  break;
367  case TGeoShape::kGeoEltu:
368  eltu = (TGeoEltu *) shape;
369  volume = 2*TMath::Pi()*eltu->GetA()*eltu->GetB()*eltu->GetDz();
370  break;
371  case TGeoShape::kGeoArb8:
372  case TGeoShape::kGeoTrap:
373  case TGeoShape::kGeoTorus:
374  case TGeoShape::kGeoPara:
375  case TGeoShape::kGeoSph:
376  case TGeoShape::kGeoCtub:
377  default:
378  if (Debug())
379  cout << bit << "\t has not yet implemented for " << shape->GetName() << endl;
380  volume = TMath::Min(volBB,volBC);
381  break;
382  }
383  break;
384  }
385  }
386  if (Debug())
387  cout << " =============> Volume = " << volume
388  << "\t" << shape->GetName() << "\tBB\t" << volBB << "\tBC\t" << volBC << endl;
389  assert(volume >= 0);
390  if (Debug() && ! (volBB - volume > -1e-7 && volBC - volume > -1.e-7) ) {
391  PrintShape(shape);
392  cout << "volume\t" << volume << "\tvolBB\t" << volBB << "\tvolBC\t" << volBC << endl;
393  cout << "\tBoundingBox dX\t" << Box->GetDX() << "\tdY\t" << Box->GetDY() << "\tdZ\t" << Box->GetDZ() << endl;
394  cout << "\tBoundingCylinder\trMin\t" << TMath::Sqrt(paramsBC[0]) << "\trMax\t" << TMath::Sqrt(paramsBC[1])
395  << "\tdZ\t" << Box->GetDZ() << endl;
396 
397  assert(volBB - volume > -1e-7 && volBC - volume > -1.e-7);
398  }
399  return volume;
400  }
401 //________________________________________________________________________________
402 Int_t StiVMCToolKit::Add2ElementList(Int_t NElem,const TGeoMaterial *mat, Elem_t *ElementList) {
403  assert(NElem>=0 && NElem<NoElemMax);
404  if (! NAir) {
405  NAir = sizeof(Air)/sizeof(ElemV_t);
406  Double_t W = 0;
407  for (Int_t i = 0; i < NAir; i++) W += Air[i].A*Air[i].V;
408  for (Int_t i = 0; i < NAir; i++) Air[i].W = Air[i].A*Air[i].V/W;
409  }
410  if (mat->InheritsFrom("TGeoMixture")) {
411  TGeoMixture *mix = (TGeoMixture *) mat;
412  Int_t N = mix->GetNelements();
413  Double_t *A = mix->GetAmixt();
414  Double_t *Z = mix->GetZmixt();
415  Double_t *W = mix->GetWmixt();
416  for (Int_t i = 0; i < N; i++) {
417  if (TMath::Abs(Z[i] - 7.30) < 1.e-3 &&
418  TMath::Abs(A[i] - 14.61) < 1.e-3) {
419  for (Int_t j = 0; j < NAir; j++) {
420  ElementList[NElem].index = NElem+1;
421  ElementList[NElem].W = W[i]*Air[j].W;
422  ElementList[NElem].A = Air[j].A;
423  ElementList[NElem].Z = Air[j].Z;
424  NElem++;
425  }
426  } else {
427  if (TMath::Abs(Z[i] - 17.436 ) < 1.e-3 &&
428  TMath::Abs(A[i] - 38.691) < 1.e-3) {
429  for (Int_t j = 0; j < NP10; j++) {
430  ElementList[NElem].index = NElem+1;
431  ElementList[NElem].W = W[i]*P10[j].W;
432  ElementList[NElem].A = P10[j].A;
433  ElementList[NElem].Z = P10[j].Z;
434  NElem++;
435  }
436  } else {
437  ElementList[NElem].index = NElem+1;
438  ElementList[NElem].W = W[i];
439  ElementList[NElem].A = A[i];
440  ElementList[NElem].Z = Z[i];
441  NElem++;
442  }
443  }
444  }
445  } else {
446  Double_t A = mat->GetA();
447  Double_t Z = mat->GetZ();
448  if (TMath::Abs(Z - 7.30) < 1.e-3 &&
449  TMath::Abs(A - 14.61) < 1.e-3) {
450  for (Int_t j = 0; j < NAir; j++) {
451  ElementList[NElem].index = NElem+1;
452  ElementList[NElem].W = Air[j].W;
453  ElementList[NElem].A = Air[j].A;
454  ElementList[NElem].Z = Air[j].Z;
455  NElem++;
456  }
457  } else {
458  if (TMath::Abs(Z - 17.436 ) < 1.e-3 &&
459  TMath::Abs(A - 38.691) < 1.e-3) {
460  for (Int_t j = 0; j < NP10; j++) {
461  ElementList[NElem].index = NElem+1;
462  ElementList[NElem].W = P10[j].W;
463  ElementList[NElem].A = P10[j].A;
464  ElementList[NElem].Z = P10[j].Z;
465  NElem++;
466  }
467  } else {
468  ElementList[NElem].index = NElem+1;
469  ElementList[NElem].W = 1.;
470  ElementList[NElem].A = A;
471  ElementList[NElem].Z = Z;
472  NElem++;
473  }
474  }
475  }
476  assert(NElem < NoElemMax);
477  return NElem;
478 }
479 //________________________________________________________________________________
480 Int_t StiVMCToolKit::Merge2ElementList(Int_t NElem, Elem_t *ElementList,
481  Int_t NElemD, Elem_t *ElementListD, Double_t weight) {
482  assert(NElem>=0 && NElem<NoElemMax);
483  for (Int_t i = 0; i < NElemD; i++) {
484  ElementList[NElem].index = NElem+1;
485  ElementList[NElem].W = weight*ElementListD[i].W;
486  ElementList[NElem].A = ElementListD[i].A;
487  ElementList[NElem].Z = ElementListD[i].Z;
488  NElem++;
489  assert(NElem < NoElemMax);
490  }
491  return NElem;
492 }
493 //________________________________________________________________________________
494 Int_t StiVMCToolKit::NormolizeElementList(Int_t NElem, Elem_t *ElementList){
495  assert(NElem>=0 && NElem<NoElemMax);
496  Double_t W = 0;
497  for (Int_t i = 0; i < NElem; i++) W += ElementList[i].W ;
498  for (Int_t i = 0; i < NElem; i++) ElementList[i].W /= W;
499  for (Int_t k = 0; k < NElem; k++) {
500  if (ElementList[k].W > 0) {
501  for (Int_t l = k+1; l < NElem; l++) {
502  if (ElementList[l].W > 0 &&
503  ElementList[k].A == ElementList[l].A &&
504  ElementList[k].Z == ElementList[l].Z) {
505  ElementList[k].W += ElementList[l].W;
506  ElementList[l].W = -1;
507  }
508  }
509  }
510  }
511  Int_t N = 0;
512 
513  for (Int_t k = 0; k < NElem; k++) {
514  if (ElementList[k].W <= 0) continue;
515  if (k != N)
516  ElementList[N] = ElementList[k];
517  N++;
518  }
519  return N;
520 }
521 //________________________________________________________________________________
522 Double_t StiVMCToolKit::GetWeight(TGeoNode *nodeT, const TString &pathT,
523  Int_t *NElem, Elem_t *ElementList) {
524  Double_t Weight = 0;
525  Double_t WeightT = 0;
526  if (! nodeT) {
527  if (! gGeoManager) return Weight;
528  gGeoManager->RestoreMasterVolume();
529  gGeoManager->CdTop();
530  gGeoManager->cd(pathT.Data());
531  nodeT = gGeoManager->GetCurrentNode();
532  if (! nodeT) return Weight;
533  }
534  TGeoVolume *volT = nodeT->GetVolume();
535  const TGeoMaterial *mat = volT->GetMaterial();
536 
537  Double_t dens = mat->GetDensity();
538  TGeoShape *shapeT = (TGeoShape *) volT->GetShape();
539  if (! shapeT) return Weight;
540  Double_t volume = GetShapeVolume(shapeT);
541  WeightT = Weight = dens*volume;
542 #if 0
543  if (! ElementList ) {
544  ElementList = new Elem_t[NoElemMax];
545  NElem = new Int_t[1]; NElem[0] = 0;
546  }
547 #else
548  assert(ElementList && NElem);
549 #endif
550  Int_t im1 = NElem[0];
551  NElem[0] = __addelement__(NElem[0], mat, ElementList);
552  Int_t im2 = NElem[0];
553  // Double_t x0 = mat->GetRadLen();
554  // cout << mat->GetName() << "\tdens = " << dens << "\tx0 = " << x0 << endl;
555  TObjArray *nodes = volT->GetNodes();
556  Int_t nd = nodeT->GetNdaughters();
557  //
558  if (nd > 0) {
559  vector<Double_t> weights(nd);
560  for (Int_t id = 0; id < nd; id++) {
561  TGeoNode *node = (TGeoNode*)nodes->UncheckedAt(id);
562  if (! node) continue;
563  // if (node->IsOverlapping()) continue;
564  TGeoVolume *vol = node->GetVolume();
565  if (! vol) continue;
566  TString name(TString(vol->GetName())); //cout << "Ask for " << name << endl;
567  // if (name.BeginsWith("FTPC")) continue;
568  TGeoShape *shapeD = (TGeoShape *) vol->GetShape();
569  if (! shapeD ) continue;
570  Double_t weight = dens*GetShapeVolume(shapeD);// cout << "\tWeight\t" << Weight << endl;
571  WeightT -= weight;
572  Weight -= weight;
573  TString path = pathT;
574  if (path != "") path += "/";
575  path += node->GetName();
576  Int_t NElemD[1] = {0};
577  Elem_t ElementListD[NoElemMax];
578  weights[id] = GetWeight(node,path,NElemD,ElementListD);
579  Weight += weights[id]; // cout << "\tWeight\t" << Weight << endl;
580  NElem[0] = Merge2ElementList(NElem[0], ElementList, NElemD[0], ElementListD, weights[id]);
581  }
582  } // nd > 0
583  for (Int_t i = im1; i < im2; i++) {// update mother weight fraction
584  ElementList[i].W *= WeightT;
585  }
586  NElem[0] = NormolizeElementList(NElem[0],ElementList);
587  return Weight;
588 }
589 
590 //________________________________________________________________________________
591 Double_t StiVMCToolKit::GetVolumeWeight(TGeoVolume *volT, Int_t *NElem, Elem_t *ElementList)
592 {
593  if (! volT || !volT->GetShape()) return 0;
594 
595  Double_t Weight = 0;
596  Double_t WeightT = 0;
597  const TGeoMaterial *mat = volT->GetMaterial();
598  TGeoShape *shapeT = (TGeoShape *) volT->GetShape();
599  if (! shapeT) return 0;
600 
601  Double_t dens = mat->GetDensity();
602  if (Debug())
603  cout << "volT \t" << volT->GetName() << "\t" << volT->GetTitle() << endl;
604  Double_t volume = GetShapeVolume(shapeT);
605  Weight = dens*volume;
606  WeightT = Weight;
607 #if 0
608  if (! ElementList ) {
609 
610  ElementList = new Elem_t[NoElemMax];
611  NElem = new Int_t[1]; NElem[0] = 0;
612  }
613 #else
614  assert(ElementList && "StiVMCToolKit::GetWeight: Create the the array first!");
615 #endif
616  Int_t im1 = *NElem;
617  Int_t im2 = __addelement__(im1, mat, ElementList);
618  *NElem = im2;
619  // Double_t x0 = mat->GetRadLen();
620  // cout << mat->GetName() << "\tdens = " << dens << "\tx0 = " << x0 << endl;
621  TObjArray *nodes = volT->GetNodes();
622  Int_t nd = volT->GetNdaughters();
623  //
624  if (nd > 0) {
625  vector<Double_t> weights(nd);
626  for (Int_t id = 0; id < nd; id++) {
627  TGeoNode *node = (TGeoNode*) nodes->UncheckedAt(id);
628  if (! node) continue;
629  // if (node->IsOverlapping()) continue;
630  TGeoVolume *vol = node->GetVolume();
631  if (! vol) continue;
632  TString name(TString(vol->GetName())); //cout << "Ask for " << name << endl;
633  // if (name.BeginsWith("FTPC")) continue;
634  TGeoShape *shapeD = (TGeoShape *) vol->GetShape();
635  if (! shapeD ) continue;
636  Double_t weight = dens*GetShapeVolume(shapeD);// cout << "\tWeight\t" << Weight << endl;
637  WeightT -= weight;
638  Weight -= weight;
639  Int_t NElemD[1] = {0};
640  Elem_t ElementListD[NoElemMax];
641  weights[id] = GetVolumeWeight(vol,NElemD,ElementListD);
642  Weight += weights[id]; // cout << "\tWeight\t" << Weight << endl;
643  Int_t im =*NElem;
644  Int_t imm = Merge2ElementList(im, ElementList, NElemD[0], ElementListD, weights[id]);
645  *NElem = imm;
646  }
647  } // nd > 0
648  for (Int_t i = im1; i < im2; i++) {// update mother weight fraction
649  ElementList[i].W *= WeightT;
650  }
651  Int_t in = *NElem;
652  Int_t inn= NormolizeElementList(in,ElementList);
653  *NElem = inn;
654  return Weight;
655 }
656 //________________________________________________________________________________
657 void StiVMCToolKit::MakeAverageVolume(TGeoVolume *volT, TGeoShape *&newshape, TGeoMedium *&newmed, Double_t *xyzM) {
658  if (! volT) return;
659  vector<Elem_t> ElementList(NoElemMax);
660  Int_t NElem = 0;
661  // NElem = ElementList.size()-1;
662  Int_t *al = &NElem;
663  Elem_t *el = &ElementList[0];
664  TGeoVolume *v = volT;
665 // Double_t Weight = GetWeight(volT, &NElem, &ElementList[0]);
666  Double_t Weight =
667  GetVolumeWeight(v
668  , al
669  , el
670  );
671  const TGeoMedium *med = volT->GetMedium();
672  const TGeoMaterial *mat = volT->GetMaterial();
673  if (Debug()) {
674  if (med->GetParam(0)) cout << "===================" << endl;
675  else cout << "+++++++++++++++++++" << endl;
676  cout << volT->GetName() << " === "
677  << " Weight " << Weight << "[g]\t";
678  cout << "material\t" << mat->GetName() << "\tmedium\t" << med->GetName() << endl;
679  }
680  TGeoShape *shapeT = (TGeoShape *) volT->GetShape();
681  Double_t volume = GetShapeVolume(shapeT);
682  if (Debug()) PrintShape(shapeT);
683  shapeT->ComputeBBox();
684  TGeoBBox * Box = (TGeoBBox *) shapeT;
685  Double_t dx=-2001, dy=-2002, dz=-2003, rmin=-2009, rmax=-2010;
686  Double_t paramsBC[4];
687  shapeT->GetBoundingCylinder(paramsBC);
688  Double_t volBB = 8*Box->GetDX()*Box->GetDY()*Box->GetDZ();
689  Double_t volBC = 2*TMath::Pi()*(paramsBC[1] - paramsBC[0])*Box->GetDZ();
690  Double_t volB = 0;
691  if (Debug()) cout << shapeT->GetName();
692  enum ShapeId {kBox = 1, kTube, kKeep};
693  Int_t iShape = 3;
694  TString name(volT->GetName());
695  if (name == "SROD") {// replace oval road by box
696  iShape = 1;
697  Double_t S = TMath::Pi()*(paramsBC[1] - paramsBC[0]);
698  dy = TMath::Sqrt(paramsBC[1]) - TMath::Sqrt(paramsBC[0]);
699  dx = S/(4*dy);
700  dz = Box->GetDZ();
701  Box->SetBoxDimensions(dx,dy,dz);
702  volB = volBB = volBC;
703  } else {
704  if ((shapeT->TestShapeBit(TGeoShape::kGeoTube) ||
705  shapeT->TestShapeBit(TGeoShape::kGeoTubeSeg)) &&
706  xyzM && xyzM[0]*xyzM[0] + xyzM[1]*xyzM[1] < 1.e-3) {
707  iShape = 3;
708  volB = volume;
709  if (Debug()) cout << "\tKeep shape ===========\t";
710  }
711  else {
712  if ( (volBB <= volBC) || (xyzM && xyzM[0]*xyzM[0] + xyzM[1]*xyzM[1] > 1.e-3) ) {
713  iShape = 1;
714  volB = volBB;
715  dx = Box->GetDX();
716  dy = Box->GetDY();
717  dz = Box->GetDZ();
718  if (Debug()) cout << "\tBoundingBox dX\t" << dx << "\tdY\t" << dy << "\tdZ\t" << dz;
719  } else {
720  iShape = 2;
721  volB = volBC;
722  rmin = TMath::Sqrt(paramsBC[0]);
723  rmax = TMath::Sqrt(paramsBC[1]);
724  dz = Box->GetDZ();
725  if (Debug()) cout << "\tBoundingCylinder\trMin\t" << rmin << "\trMax\t" << rmax << "\tdZ\t" << dz;
726  }
727  }
728  }
729  Double_t Dens = Weight/volB;
730  if (Debug()) cout << "\tvolume\t" << volume << "\tvolB\t" << volB << "\tDens\t" << Dens << endl;
731  TString newmatname(mat->GetName());
732  newmatname += "_";
733  newmatname += name;
734  TGeoMixture *newmat = new TGeoMixture(newmatname.Data(),NElem, Dens);
735  for (Int_t i = 0; i < NElem; i++) {
736  if (Debug())
737  cout << Form("id%3i\tW%10.3g\tA%5.1f\t%5.1f\n",
738  ElementList[i].index,ElementList[i].W,ElementList[i].A,ElementList[i].Z);
739  // cout << "id\t" << ElementList[i].index << "\tW\t" << ElementList[i].W
740  // << "\tA\t" << ElementList[i].A << "\tZ\t" << ElementList[i].Z << endl;
741  newmat->DefineElement(i,ElementList[i].A,ElementList[i].Z,ElementList[i].W);
742  }
743  newmat->SetUniqueID(gGeoManager->AddMaterial(newmat));
744  Int_t matid = newmat->GetUniqueID();
745  // name numed imat,
746  Double_t params[10];
747  for (Int_t i = 0; i < 10; i++) params[i] = med->GetParam(i);
748  newmed = new TGeoMedium(newmatname.Data(),matid,newmat,params);
749  if (Debug() && xyzM) {
750  cout << "\txyzM x\t" << xyzM[0] << "\ty\t" << xyzM[1] << "\tz\t" << xyzM[2] << endl;
751  }
752  newshape = shapeT;
753  if (iShape == 1) newshape = new TGeoBBox(dx, dy, dz);
754  else if (iShape == 2) {
755  newshape = new TGeoTube(rmin, rmax, dz);
756  }
757 }
758 //________________________________________________________________________________
759 TGeoManager * StiVMCToolKit::GetVMC() {
760  TGeoManager *gGeo = 0;
761 #ifndef __NOVMC__
762 
764  gGeo = gGeoManager;
765  if (gGeo) return gGeo;
766  LOG_INFO << "StiVMCToolKit::GetVMC() : Get VMC geometry" <<endm;
767  if (StMaker::GetChain()) {
768  StMaker::GetChain()->GetDataBase("VmcGeometry");
769  }
770  if (! gGeoManager)
771  LOG_ERROR << "StiVMCToolKit::GetVMC() : Can't get VMC geometry" <<endm;
772  gGeo = gGeoManager;
773 #endif
774  return gGeo;
775 }
776 //________________________________________________________________________________
777 void StiVMCToolKit::TestVMC4Reconstruction(){
780  GetVMC();
781  // TObjArray *volumes = gGeoManager->GetListOfVolumes();
782  Int_t NV[2] = {nTopVol,NofVolToBEAveraged};
783  VolumeMap_t *list = 0;
784  TGeoVolume *volT = 0;
785  cout << "<table>" << endl;
786  cout << "<tr><td>name</td><td>Comment </td><td>Volume[cm**3]</td><td>Weight[g]</td><td>Av.Dens.[g/cm**3]</td></tr>" << endl;
787  for (Int_t i = 0; i < 2; i++) {
788  if (i == 0) list = TopVolumes;
789  else list = VolumesToBeAveraged;
790  for (Int_t j = 0; j < NV[i]; j++) {
791  volT = gGeoManager->GetVolume(list[j].name);
792  // volT = (TGeoVolume *) volumes->FindObject(list[j].name);
793  if (! volT) {cout << "Can't find " << list[j].name << "\t" << list[j].comment << endl; continue;}
794  TGeoShape *shapeT = volT->GetShape();
795  Double_t volume = GetShapeVolume(shapeT);
796  Double_t weight = GetVolumeWeight(volT, 0, 0); // it leads to the memory leak !!!
797  cout << "<tr><td>"<< list[j].name << "</td><td>" << list[j].comment << "</td>"
798  << "<td>" << Form("%10.3g",volume) << "</td>"
799  << "<td>" << Form("%10.3g",weight) << "</td>"
800  << "<td>" << Form("%10.3g",weight/volume) << "</td></tr>"
801  << endl;
802  }
803  }
804  cout << "</table>" << endl;
805 }
806 //________________________________________________________________________________
807 TGeoPhysicalNode *StiVMCToolKit::Alignment(const TGeoNode *nodeT, const Char_t *pathT,
808  TGeoVolume *volT, TGeoShape *newshape, TGeoMedium* newmed) {
809  // Alignment
810  TObjArray *listP = gGeoManager->GetListOfPhysicalNodes();
811  TGeoPhysicalNode *nodeP = 0;
812  TGeoCombiTrans *trP = 0;
813  if (listP) {
814  Int_t N = listP->GetEntries();
815  for (Int_t i = 0; i < N; i++) {
816  TGeoPhysicalNode *nod = dynamic_cast<TGeoPhysicalNode *> (listP->At(i));
817  if (nod && TString(nod->GetName()) == TString(pathT)) {
818  nodeP = nod; break;
819  }
820  }
821  }
822  Double_t local[3] = {0,0,0};
823  TGeoShape *shapeT = volT->GetShape();
824  if (Debug()) {
825  LOG_INFO << "StiVMCToolKit::Alignment : node\t" << nodeT->GetName()
826  << "\tvolume\t" << volT->GetName() << "\t:" << volT->GetTitle()
827  << "\tshape\t" << shapeT->GetName() << "\tmaterial\t" << volT->GetMaterial()->GetName();
828  if (newshape) LOG_INFO << "\tnewshape\t" << newshape->GetName();
829  LOG_INFO << endm;
830  }
831  if (shapeT->TestShapeBit(TGeoShape::kGeoPcon) || shapeT->TestShapeBit(TGeoShape::kGeoPgon)) {
832  TGeoPcon *pcon = (TGeoPcon *) shapeT;
833  Int_t Nz = pcon->GetNz();
834  local[2] = 0.5*(pcon->GetZ(0) + pcon->GetZ(Nz-1));
835  }
836  Double_t master[3];
837  nodeT->LocalToMaster(local,master);
838  if (Debug())
839  cout << "\tmaster x\t" << master[0] << "\ty\t" << master[1] << "\tz\t" << master[2] << endl;
840  if (!nodeP) nodeP = gGeoManager->MakePhysicalNode(pathT);
841  if (!nodeP) return nodeP;
842 #if 0
843  if (nodeP->IsAligned()) {
844  trP = nodeP->GetNode()->GetMatrix();
845  trP->SetTranslation(master[0], master[1], master[2]);
846  } else {
847  trP = new TGeoCombiTrans(master[0], master[1], master[2], nodeP->GetNode()->GetMatrix());
848  }
849 #endif
850  nodeP->Align(trP,newshape);//,kTRUE);
851  TGeoVolume *newvol = nodeP->GetNode(-1)->GetVolume();
852  if (Debug())
853  cout << "newvol\t" << newvol->GetName() << "\tmed\t" << newvol->GetMedium()->GetName() << endl;
854  newvol->SetMedium(newmed);
855  newvol->SetTitle("AVERAGED");
856  newvol->SetLineColor(1);
857  newvol->SetVisibility(1);
858  TObjArray *nodes = newvol->GetNodes();
859  if (nodes && ! newvol->TObject::TestBit(TGeoVolume::kVolumeImportNodes)) {
860  nodes->Delete();
861  delete nodes;
862  newvol->ClearNodes();
863  }
864  return nodeP;
865 }
866 //________________________________________________________________________________
867 TGeoPhysicalNode *StiVMCToolKit::LoopOverNodes(const TGeoNode *nodeT, const Char_t *pathT,
868  const Char_t *name,void ( *callback)(TGeoPhysicalNode *)){
869  TGeoVolume *volT = nodeT->GetVolume();
870  const Char_t *nameT = volT->GetName();
871  TGeoPhysicalNode *nodeP = 0;
872  TGeoShape *newshape = 0;
873  TGeoMedium* newmed = 0;
874  if (nameT && name && ! strncmp(nameT,name,4)) {
875  // if (TString(nameT) == TString(VolumesToBeAveraged[i].name)) {
876  Double_t local[3] = {0, 0, 0};
877  TGeoShape *shapeT = volT->GetShape();
878  if (shapeT->TestShapeBit(TGeoShape::kGeoPcon) || shapeT->TestShapeBit(TGeoShape::kGeoPgon)) {
879  TGeoPcon *pcon = (TGeoPcon *) shapeT;
880  Int_t Nz = pcon->GetNz();
881  local[2] = 0.5*(pcon->GetZ(0) + pcon->GetZ(Nz-1));
882  }
883  Double_t master[3];
884  nodeT->LocalToMaster(local,master);
885  MakeAverageVolume(volT, newshape, newmed, master);
886  nodeP = Alignment(nodeT,pathT, volT, newshape, newmed);
887  if (nodeP) {
888  if (! callback) MakeVolume(nodeP);
889  else callback(nodeP);
890  }
891  return nodeP;
892  }
893  TObjArray *nodes = volT->GetNodes();
894  Int_t nd = volT->GetNdaughters();
895  for (Int_t id = 0; id < nd; id++) {
896  TGeoNode *node = (TGeoNode*) nodes->UncheckedAt(id);
897  if (! node) continue;
898  TString path = pathT;
899  if (path != "") path += "/";
900  path += node->GetName();
901  if (! name && Debug()) {
902  cout << path;
903  TGeoVolume *vol = node->GetVolume();
904  const TGeoMedium *med = vol->GetMedium();
905  if (med->GetParam(0)) {cout << "\t===================" << endl; continue;}
906  else cout << endl;
907  }
908  LoopOverNodes(node, path, name, callback);
909  }
910  return nodeP;
911 }
912 //________________________________________________________________________________
913 void StiVMCToolKit::MakeVolume(TGeoPhysicalNode *nodeP) {
914  if (Debug()) {
915  LOG_INFO << "StiVMCToolKit::MakeVolume : TGeoPhysicalNode\t" << nodeP->GetName() << endm;
916  TGeoVolume *volP = nodeP->GetVolume();
917  TGeoMaterial *matP = volP->GetMaterial(); matP->Print("");
918  TGeoShape *shapeP = nodeP->GetShape(); cout << "New Shape\t"; PrintShape(shapeP);
919  TGeoHMatrix *hmat = nodeP->GetMatrix(); hmat->Print("");
920  }
921 }
922 //________________________________________________________________________________
923 Double_t StiVMCToolKit::GetPotI(const TGeoMaterial *mat) {
924  Double_t PotI = 0;
925  if (mat) {
926  Double_t s1 = 0, s2 = 0;
927  if (mat->InheritsFrom("TGeoMixture")) {
928  TGeoMixture *mix = (TGeoMixture *) mat;
929  Int_t N = mix->GetNelements();
930  assert(N);
931  Double_t *A = mix->GetAmixt();
932  Double_t *Z = mix->GetZmixt();
933  Double_t *W = mix->GetWmixt();
934  for (Int_t i = 0; i < N; i++) {
935  s1 += W[i]*Z[i]/A[i];
936  s2 += W[i]*Z[i]*TMath::Log(Z[i])/A[i];
937  }
938  PotI=16.e-9*TMath::Exp(0.9*s2/s1);
939  } else {
940  Double_t Z = mat->GetZ();
941  PotI=16.e-9*TMath::Power(Z,0.9);
942  }
943  }
944  return PotI;
945 }
946 //________________________________________________________________________________
947 void StiVMCToolKit::GetVMC4Reconstruction(const Char_t *pathT, const Char_t *nameT){
953  GetVMC();
954  gGeoManager->RestoreMasterVolume();
955  gGeoManager->CdTop();
956  TString path("");
957  if (pathT) {gGeoManager->cd(pathT); path = pathT;}
958  else {path = gGeoManager->GetCurrentNode()->GetName();}
959  TGeoNode *nodeT = gGeoManager->GetCurrentNode();
960  if (! nodeT) return;
961 #if 1
962  LoopOverNodes(nodeT, path, nameT);
963 #else
964  for (Int_t i = 0; i < NofVolToBEAveraged; i++) {
965  LoopOverNodes(nodeT, path, VolumesToBeAveraged[i].name);
966  }
967 #endif
968 }
969 #ifndef __ROOT__
970 //________________________________________________________________________________
971 void TestVMCTK() {
974  // SVT weight
975  TString path("HALL_1/CAVE_1/SVTT_1");
976  Double_t weight = 1.e-3*StiVMCToolKit::GetWeight(0,path, 0, 0); // It leads to the memory leak
977  cout << "StiVMCToolKit::TestVMCTK() -I- total weight for "
978  << path.Data() << "\t" << weight << "[kg]" << endl;
979  StiVMCToolKit::TestVMC4Reconstruction();
980 }
981 #endif
982 //________________________________________________________________________________
983 Double_t StiVMCToolKit::Nice(Double_t phi) {
984  while (phi < 2*TMath::Pi()) phi += 2*TMath::Pi();
985  while (phi >= 2*TMath::Pi()) phi -= 2*TMath::Pi();
986  return phi;
987 }
VolumeMap_t
Definition: StiVMCToolKit.h:28
Elem_t
Definition: StiVMCToolKit.cxx:22
ElemV_t
Definition: StiVMCToolKit.cxx:29
Generated by   doxygen 1.8.5