]>
Commit | Line | Data |
---|---|---|
e2e989f9 | 1 | /************************************************************************** |
2 | * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * | |
3 | * * | |
4 | * Author: The ALICE Off-line Project. * | |
5 | * Contributors are mentioned in the code where appropriate. * | |
6 | * * | |
7 | * Permission to use, copy, modify and distribute this software and its * | |
8 | * documentation strictly for non-commercial purposes is hereby granted * | |
9 | * without fee, provided that the above copyright notice appears in all * | |
10 | * copies and that both the copyright notice and this permission notice * | |
11 | * appear in the supporting documentation. The authors make no claims * | |
12 | * about the suitability of this software for any purpose. It is * | |
13 | * provided "as is" without express or implied warranty. * | |
14 | **************************************************************************/ | |
15 | ||
16 | // $Id$ | |
c1a0a2f0 | 17 | |
18 | // Class TFlukaMCGeometry | |
19 | // -------------------- | |
20 | // Implementation of the TVirtualMCGeometry interface | |
21 | // for defining and using TGeo geometry with FLUKA. | |
22 | // This allows the FLUKA MonteCarlo to run with the TGeo | |
23 | // geometrical modeller | |
8495a208 | 24 | // Author: Andrei Gheata 10/07/2003 |
25 | ||
c1a0a2f0 | 26 | #include "Riostream.h" |
27 | ||
28 | #include "TCallf77.h" | |
89125dc7 | 29 | #include "TFluka.h" |
8495a208 | 30 | #include "TFlukaMCGeometry.h" |
31 | #include "TGeoManager.h" | |
32 | #include "TGeoVolume.h" | |
c1a0a2f0 | 33 | #include "TObjString.h" |
8495a208 | 34 | |
35 | #ifndef WIN32 | |
36 | # define idnrwr idnrwr_ | |
37 | # define g1wr g1wr_ | |
38 | # define g1rtwr g1rtwr_ | |
39 | # define conhwr conhwr_ | |
40 | # define inihwr inihwr_ | |
41 | # define jomiwr jomiwr_ | |
42 | # define lkdbwr lkdbwr_ | |
43 | # define lkfxwr lkfxwr_ | |
44 | # define lkmgwr lkmgwr_ | |
45 | # define lkwr lkwr_ | |
46 | # define magfld magfld_ | |
47 | # define nrmlwr nrmlwr_ | |
48 | # define rgrpwr rgrpwr_ | |
49 | # define isvhwr isvhwr_ | |
50 | ||
51 | #else | |
52 | ||
53 | # define idnrwr IDNRWR | |
54 | # define g1wr G1WR | |
55 | # define g1rtwr G1RTWR | |
56 | # define conhwr CONHWR | |
57 | # define inihwr INIHWR | |
58 | # define jomiwr JOMIWR | |
59 | # define lkdbwr LKDBWR | |
60 | # define lkfxwr LKFXWR | |
61 | # define lkmgwr LKMGWR | |
62 | # define lkwr LKWR | |
63 | # define magfld MAGFLD | |
64 | # define nrmlwr NRMLWR | |
65 | # define rgrpwr RGRPWR | |
66 | # define isvhwr ISVHWR | |
67 | ||
68 | #endif | |
69 | ||
70 | //____________________________________________________________________________ | |
71 | extern "C" | |
72 | { | |
73 | // | |
74 | // Prototypes for FLUKA navigation methods | |
75 | // | |
76 | Int_t type_of_call idnrwr(const Int_t & /*nreg*/, const Int_t & /*mlat*/); | |
77 | void type_of_call g1wr(Double_t & /*pSx*/, Double_t & /*pSy*/, Double_t & /*pSz*/, | |
78 | Double_t * /*pV*/, Int_t & /*oldReg*/ , const Int_t & /*oldLttc*/, Double_t & /*propStep*/, | |
79 | Int_t & /*nascFlag*/, Double_t & /*retStep*/, Int_t & /*newReg*/, | |
80 | Double_t & /*saf*/, Int_t & /*newLttc*/, Int_t & /*LttcFlag*/, | |
81 | Double_t *s /*Lt*/, Int_t * /*jrLt*/); | |
82 | ||
83 | void type_of_call g1rtwr(); | |
84 | void type_of_call conhwr(Int_t & /*intHist*/, Int_t * /*incrCount*/); | |
85 | void type_of_call inihwr(Int_t & /*intHist*/); | |
86 | void type_of_call jomiwr(const Int_t & /*nge*/, const Int_t & /*lin*/, const Int_t & /*lou*/, | |
87 | Int_t & /*flukaReg*/); | |
88 | void type_of_call lkdbwr(Double_t & /*pSx*/, Double_t & /*pSy*/, Double_t & /*pSz*/, | |
89 | Double_t * /*pV*/, const Int_t & /*oldReg*/, const Int_t & /*oldLttc*/, | |
90 | Int_t & /*newReg*/, Int_t & /*flagErr*/, Int_t & /*newLttc*/); | |
91 | void type_of_call lkfxwr(Double_t & /*pSx*/, Double_t & /*pSy*/, Double_t & /*pSz*/, | |
92 | Double_t * /*pV*/, const Int_t & /*oldReg*/, const Int_t & /*oldLttc*/, | |
93 | Int_t & /*newReg*/, Int_t & /*flagErr*/, Int_t & /*newLttc*/); | |
94 | void type_of_call lkmgwr(Double_t & /*pSx*/, Double_t & /*pSy*/, Double_t & /*pSz*/, | |
95 | Double_t * /*pV*/, const Int_t & /*oldReg*/, const Int_t & /*oldLttc*/, | |
96 | Int_t & /*flagErr*/, Int_t & /*newReg*/, Int_t & /*newLttc*/); | |
97 | void type_of_call lkwr(Double_t & /*pSx*/, Double_t & /*pSy*/, Double_t & /*pSz*/, | |
98 | Double_t * /*pV*/, const Int_t & /*oldReg*/, const Int_t & /*oldLttc*/, | |
99 | Int_t & /*newReg*/, Int_t & /*flagErr*/, Int_t & /*newLttc*/); | |
efde9b4d | 100 | // void type_of_call magfld(const Double_t & /*pX*/, const Double_t & /*pY*/, const Double_t & /*pZ*/, |
101 | // Double_t & /*cosBx*/, Double_t & /*cosBy*/, Double_t & /*cosBz*/, | |
102 | // Double_t & /*Bmag*/, Int_t & /*reg*/, Int_t & /*idiscflag*/); | |
8495a208 | 103 | void type_of_call nrmlwr(Double_t & /*pSx*/, Double_t & /*pSy*/, Double_t & /*pSz*/, |
104 | Double_t & /*pVx*/, Double_t & /*pVy*/, Double_t & /*pVz*/, | |
105 | Double_t * /*norml*/, const Int_t & /*oldReg*/, | |
106 | const Int_t & /*newReg*/, Int_t & /*flagErr*/); | |
107 | void type_of_call rgrpwr(const Int_t & /*flukaReg*/, const Int_t & /*ptrLttc*/, Int_t & /*g4Reg*/, | |
108 | Int_t * /*indMother*/, Int_t * /*repMother*/, Int_t & /*depthFluka*/); | |
109 | Int_t type_of_call isvhwr(const Int_t & /*fCheck*/, const Int_t & /*intHist*/); | |
110 | }; | |
111 | ||
112 | // TFluka global pointer | |
89125dc7 | 113 | TFluka *gFluka = 0; |
114 | TFlukaMCGeometry *gMCGeom = 0; | |
6f5667d1 | 115 | Int_t gNstep = 0; |
8495a208 | 116 | |
117 | ClassImp(TFlukaMCGeometry) | |
118 | ||
119 | TFlukaMCGeometry* TFlukaMCGeometry::fgInstance=0; | |
120 | ||
121 | //_____________________________________________________________________________ | |
122 | TFlukaMCGeometry::TFlukaMCGeometry(const char *name, const char *title) | |
89125dc7 | 123 | : TNamed(name, title) |
8495a208 | 124 | { |
125 | // | |
126 | // Standard constructor | |
127 | // | |
2bc4c610 | 128 | fDebug = kFALSE; |
efde9b4d | 129 | fLastMaterial = 0; |
b1536e91 | 130 | fNextRegion = 0; |
131 | fNextLattice = 0; | |
132 | fRegionList = 0; | |
89125dc7 | 133 | gFluka = (TFluka*)gMC; |
134 | gMCGeom = this; | |
6f5667d1 | 135 | gNstep = 0; |
89125dc7 | 136 | fMatList = new TObjArray(256); |
137 | fMatNames = new TObjArray(256); | |
8495a208 | 138 | } |
139 | ||
140 | //_____________________________________________________________________________ | |
141 | TFlukaMCGeometry::TFlukaMCGeometry() | |
8495a208 | 142 | { |
143 | // | |
144 | // Default constructor | |
145 | // | |
2bc4c610 | 146 | fDebug = kFALSE; |
efde9b4d | 147 | fLastMaterial = 0; |
b1536e91 | 148 | fNextRegion = 0; |
149 | fNextLattice = 0; | |
150 | fRegionList = 0; | |
89125dc7 | 151 | gFluka = (TFluka*)gMC; |
152 | gMCGeom = this; | |
6f5667d1 | 153 | gNstep = 0; |
89125dc7 | 154 | fMatList = 0; |
155 | fMatNames = 0; | |
8495a208 | 156 | } |
157 | ||
158 | //_____________________________________________________________________________ | |
159 | TFlukaMCGeometry::~TFlukaMCGeometry() | |
160 | { | |
161 | // | |
162 | // Destructor | |
163 | // | |
164 | fgInstance=0; | |
b1536e91 | 165 | if (fRegionList) delete [] fRegionList; |
89125dc7 | 166 | if (fMatList) delete fMatList; |
167 | if (fMatNames) {fMatNames->Delete(); delete fMatNames;} | |
8495a208 | 168 | if (gGeoManager) delete gGeoManager; |
169 | } | |
170 | ||
171 | // | |
172 | // private methods | |
173 | // | |
174 | //_____________________________________________________________________________ | |
175 | TFlukaMCGeometry::TFlukaMCGeometry(const TFlukaMCGeometry &) | |
89125dc7 | 176 | : TNamed() |
8495a208 | 177 | { |
178 | // | |
179 | // Copy constructor | |
180 | // | |
181 | } | |
182 | ||
183 | //_____________________________________________________________________________ | |
184 | Double_t* TFlukaMCGeometry::CreateDoubleArray(Float_t* array, Int_t size) const | |
185 | { | |
186 | // Converts Float_t* array to Double_t*, | |
187 | // !! The new array has to be deleted by user. | |
188 | // --- | |
189 | ||
190 | Double_t* doubleArray; | |
191 | if (size>0) { | |
192 | doubleArray = new Double_t[size]; | |
193 | for (Int_t i=0; i<size; i++) doubleArray[i] = array[i]; | |
194 | } | |
195 | else { | |
196 | //doubleArray = 0; | |
197 | doubleArray = new Double_t[1]; | |
198 | } | |
199 | return doubleArray; | |
200 | } | |
201 | // | |
202 | // public methods | |
8495a208 | 203 | |
8495a208 | 204 | |
8495a208 | 205 | //_____________________________________________________________________________ |
206 | Int_t TFlukaMCGeometry::GetMedium() const | |
207 | { | |
208 | // Get current medium number | |
209 | Int_t imed = 0; | |
210 | TGeoNode *node = gGeoManager->GetCurrentNode(); | |
211 | if (!node) imed = gGeoManager->GetTopNode()->GetVolume()->GetMedium()->GetId(); | |
212 | else imed = node->GetVolume()->GetMedium()->GetId(); | |
2bc4c610 | 213 | if (fDebug) printf("GetMedium=%i\n", imed); |
8495a208 | 214 | return imed; |
215 | } | |
216 | ||
0bc73b6c | 217 | //_____________________________________________________________________________ |
218 | Int_t TFlukaMCGeometry::GetFlukaMaterial(Int_t imed) const | |
219 | { | |
220 | // Returns FLUKA material index for medium IMED | |
221 | TGeoMedium *med = (TGeoMedium*)gGeoManager->GetListOfMedia()->At(imed-1); | |
222 | if (!med) { | |
223 | Error("GetFlukaMaterial", "MEDIUM %i nor found", imed); | |
224 | return -1; | |
225 | } | |
226 | Int_t imatfl = med->GetMaterial()->GetIndex(); | |
227 | return imatfl; | |
228 | } | |
229 | ||
b1536e91 | 230 | //_____________________________________________________________________________ |
231 | Int_t *TFlukaMCGeometry::GetRegionList(Int_t imed, Int_t &nreg) | |
232 | { | |
233 | // Get an ordered list of regions matching a given medium number | |
234 | nreg = 0; | |
235 | if (!fRegionList) fRegionList = new Int_t[NofVolumes()+1]; | |
236 | TIter next(gGeoManager->GetListOfUVolumes()); | |
237 | TGeoVolume *vol; | |
238 | Int_t imedium, ireg; | |
239 | while ((vol = (TGeoVolume*)next())) { | |
240 | imedium = vol->GetMedium()->GetId(); | |
241 | if (imedium == imed) { | |
242 | ireg = vol->GetNumber(); | |
243 | fRegionList[nreg++] = ireg; | |
244 | } | |
245 | } | |
246 | return fRegionList; | |
247 | } | |
248 | ||
249 | //_____________________________________________________________________________ | |
250 | Int_t *TFlukaMCGeometry::GetMaterialList(Int_t imat, Int_t &nreg) | |
251 | { | |
252 | // Get an ordered list of regions matching a given medium number | |
253 | nreg = 0; | |
254 | if (!fRegionList) fRegionList = new Int_t[NofVolumes()+1]; | |
255 | TIter next(gGeoManager->GetListOfUVolumes()); | |
256 | TGeoVolume *vol; | |
257 | Int_t imaterial, ireg; | |
258 | while ((vol = (TGeoVolume*)next())) { | |
259 | imaterial = vol->GetMedium()->GetMaterial()->GetIndex(); | |
260 | if (imaterial == imat) { | |
261 | ireg = vol->GetNumber(); | |
262 | fRegionList[nreg++] = ireg; | |
263 | } | |
264 | } | |
265 | return fRegionList; | |
266 | } | |
8495a208 | 267 | |
8495a208 | 268 | //_____________________________________________________________________________ |
269 | Int_t TFlukaMCGeometry::NofVolumes() const | |
270 | { | |
271 | // | |
272 | // Return total number of volumes in the geometry | |
273 | // | |
274 | ||
275 | return gGeoManager->GetListOfUVolumes()->GetEntriesFast()-1; | |
276 | } | |
8495a208 | 277 | |
278 | //_____________________________________________________________________________ | |
89125dc7 | 279 | TGeoMaterial * TFlukaMCGeometry::GetMakeWrongMaterial(Double_t z) |
280 | { | |
281 | // Try to replace a wrongly-defined material | |
282 | static Double_t kz[23] = {7.3, 17.8184, 7.2167, 10.856, 8.875, 8.9, 7.177, | |
283 | 25.72, 6.2363, 7.1315, 47.7056, 10.6467, 7.8598, 2.10853, 10.6001, 9.1193, | |
284 | 15.3383, 4.55, 9.6502, 6.4561, 21.7963, 29.8246, 15.4021}; | |
285 | ||
286 | Int_t ind; | |
287 | Double_t dz; | |
288 | for (ind=0; ind<23; ind++) { | |
289 | dz = TMath::Abs(z-kz[ind]); | |
290 | if (dz<1E-4) break; | |
291 | } | |
292 | if (ind>22) { | |
293 | printf("Cannot patch material with Z=%g\n", z); | |
294 | return 0; | |
295 | } | |
296 | TGeoMixture *mix = 0; | |
297 | TGeoElement *element; | |
298 | TGeoElementTable *table = TGeoElementTable::Instance(); | |
299 | switch (ind) { | |
300 | case 0: // AIR | |
301 | mix = new TGeoMixture("AIR", 4, 0.001205); | |
302 | element = table->GetElement(6); // C | |
303 | mix->DefineElement(0, element, 0.000124); | |
304 | element = table->GetElement(7); // N | |
305 | mix->DefineElement(1, element, 0.755267); | |
306 | element = table->GetElement(8); // O | |
307 | mix->DefineElement(2, element, 0.231781); | |
308 | element = table->GetElement(18); // AR | |
309 | mix->DefineElement(3, element, 0.012827); | |
310 | break; | |
311 | case 1: //SDD SI CHIP | |
312 | mix = new TGeoMixture("SDD_SI", 6, 2.4485); | |
313 | element = table->GetElement(1); | |
314 | mix->DefineElement(0, element, 0.004367771); | |
315 | element = table->GetElement(6); | |
316 | mix->DefineElement(1, element, 0.039730642); | |
317 | element = table->GetElement(7); | |
318 | mix->DefineElement(2, element, 0.001396798); | |
319 | element = table->GetElement(8); | |
320 | mix->DefineElement(3, element, 0.01169634); | |
321 | element = table->GetElement(14); | |
322 | mix->DefineElement(4, element, 0.844665); | |
323 | element = table->GetElement(47); | |
324 | mix->DefineElement(5, element, 0.09814344903); | |
325 | break; | |
326 | case 2: // WATER | |
327 | mix = new TGeoMixture("WATER", 2, 1.0); | |
328 | element = table->GetElement(1); | |
329 | mix->DefineElement(0, element, 0.111898344); | |
330 | element = table->GetElement(8); | |
331 | mix->DefineElement(1, element, 0.888101656); | |
332 | break; | |
333 | case 3: // CERAMICS | |
334 | mix = new TGeoMixture("CERAMICS", 5, 3.6); | |
335 | element = table->GetElement(8); | |
336 | mix->DefineElement(0, element, 0.59956); | |
337 | element = table->GetElement(13); | |
338 | mix->DefineElement(1, element, 0.3776); | |
339 | element = table->GetElement(14); | |
340 | mix->DefineElement(2, element, 0.00933); | |
341 | element = table->GetElement(24); | |
342 | mix->DefineElement(3, element, 0.002); | |
343 | element = table->GetElement(25); | |
344 | mix->DefineElement(4, element, 0.0115); | |
345 | break; | |
346 | case 4: // EPOXY | |
347 | mix = new TGeoMixture("G10FR4", 4, 1.8); | |
348 | element = table->GetElement(1); | |
349 | mix->DefineElement(0, element, 0.19); | |
350 | element = table->GetElement(6); | |
351 | mix->DefineElement(1, element, 0.18); | |
352 | element = table->GetElement(8); | |
353 | mix->DefineElement(2, element, 0.35); | |
354 | element = table->GetElement(14); | |
355 | mix->DefineElement(3, element, 0.28); | |
356 | break; | |
357 | case 5: // EPOXY | |
358 | mix = new TGeoMixture("G10FR4", 4, 1.8); | |
359 | element = table->GetElement(1); | |
360 | mix->DefineElement(0, element, 0.19); | |
361 | element = table->GetElement(6); | |
362 | mix->DefineElement(1, element, 0.18); | |
363 | element = table->GetElement(8); | |
364 | mix->DefineElement(2, element, 0.35); | |
365 | element = table->GetElement(14); | |
366 | mix->DefineElement(3, element, 0.28); | |
367 | break; | |
368 | case 6: // KAPTON | |
369 | mix = new TGeoMixture("KAPTON", 4, 1.3); | |
370 | element = table->GetElement(1); | |
371 | mix->DefineElement(0, element, 0.026363415); | |
372 | element = table->GetElement(6); | |
373 | mix->DefineElement(1, element, 0.6911272); | |
374 | element = table->GetElement(7); | |
375 | mix->DefineElement(2, element, 0.073271325); | |
376 | element = table->GetElement(8); | |
377 | mix->DefineElement(3, element, 0.209238060); | |
378 | break; | |
379 | case 7: // INOX | |
380 | mix = new TGeoMixture("INOX", 9, 7.9); | |
381 | element = table->GetElement(6); | |
382 | mix->DefineElement(0, element, 0.0003); | |
383 | element = table->GetElement(14); | |
384 | mix->DefineElement(1, element, 0.01); | |
385 | element = table->GetElement(15); | |
386 | mix->DefineElement(2, element, 0.00045); | |
387 | element = table->GetElement(16); | |
388 | mix->DefineElement(3, element, 0.0003); | |
389 | element = table->GetElement(24); | |
390 | mix->DefineElement(4, element, 0.17); | |
391 | element = table->GetElement(25); | |
392 | mix->DefineElement(5, element, 0.02); | |
393 | element = table->GetElement(26); | |
394 | mix->DefineElement(6, element, 0.654); | |
395 | element = table->GetElement(28); | |
396 | mix->DefineElement(7, element, 0.12); | |
397 | element = table->GetElement(42); | |
398 | mix->DefineElement(8, element, 0.025); | |
399 | break; | |
400 | case 8: // ROHACELL | |
401 | mix = new TGeoMixture("ROHACELL", 4, 0.05); | |
402 | element = table->GetElement(1); | |
403 | mix->DefineElement(0, element, 0.07836617); | |
404 | element = table->GetElement(6); | |
405 | mix->DefineElement(1, element, 0.64648941); | |
406 | element = table->GetElement(7); | |
407 | mix->DefineElement(2, element, 0.08376983); | |
408 | element = table->GetElement(8); | |
409 | mix->DefineElement(3, element, 0.19137459); | |
410 | break; | |
411 | case 9: // SDD-C-AL | |
412 | mix = new TGeoMixture("SDD-C-AL", 5, 1.9837); | |
413 | element = table->GetElement(1); | |
414 | mix->DefineElement(0, element, 0.022632); | |
415 | element = table->GetElement(6); | |
416 | mix->DefineElement(1, element, 0.8176579); | |
417 | element = table->GetElement(7); | |
418 | mix->DefineElement(2, element, 0.0093488); | |
419 | element = table->GetElement(8); | |
420 | mix->DefineElement(3, element, 0.0503618); | |
421 | element = table->GetElement(13); | |
422 | mix->DefineElement(4, element, 0.1); | |
423 | break; | |
424 | case 10: // X7R-CAP | |
425 | mix = new TGeoMixture("X7R-CAP", 7, 6.72); | |
426 | element = table->GetElement(8); | |
427 | mix->DefineElement(0, element, 0.085975822); | |
428 | element = table->GetElement(22); | |
429 | mix->DefineElement(1, element, 0.084755042); | |
430 | element = table->GetElement(28); | |
431 | mix->DefineElement(2, element, 0.038244751); | |
432 | element = table->GetElement(29); | |
433 | mix->DefineElement(3, element, 0.009471271); | |
434 | element = table->GetElement(50); | |
435 | mix->DefineElement(4, element, 0.321736471); | |
436 | element = table->GetElement(56); | |
437 | mix->DefineElement(5, element, 0.251639432); | |
438 | element = table->GetElement(82); | |
439 | mix->DefineElement(6, element, 0.2081768); | |
440 | break; | |
441 | case 11: // SDD ruby sph. Al2O3 | |
442 | mix = new TGeoMixture("AL2O3", 2, 3.97); | |
443 | element = table->GetElement(8); | |
444 | mix->DefineElement(0, element, 0.5293); | |
445 | element = table->GetElement(13); | |
446 | mix->DefineElement(1, element, 0.4707); | |
447 | break; | |
448 | case 12: // SDD HV microcable | |
449 | mix = new TGeoMixture("HV-CABLE", 5, 1.6087); | |
450 | element = table->GetElement(1); | |
451 | mix->DefineElement(0, element, 0.01983871336); | |
452 | element = table->GetElement(6); | |
453 | mix->DefineElement(1, element, 0.520088819984); | |
454 | element = table->GetElement(7); | |
455 | mix->DefineElement(2, element, 0.0551367996); | |
456 | element = table->GetElement(8); | |
457 | mix->DefineElement(3, element, 0.157399667056); | |
458 | element = table->GetElement(13); | |
459 | mix->DefineElement(4, element, 0.247536); | |
460 | break; | |
461 | case 13: //SDD LV+signal cable | |
462 | mix = new TGeoMixture("LV-CABLE", 5, 2.1035); | |
463 | element = table->GetElement(1); | |
464 | mix->DefineElement(0, element, 0.0082859922); | |
465 | element = table->GetElement(6); | |
466 | mix->DefineElement(1, element, 0.21722436468); | |
467 | element = table->GetElement(7); | |
468 | mix->DefineElement(2, element, 0.023028867); | |
469 | element = table->GetElement(8); | |
470 | mix->DefineElement(3, element, 0.06574077612); | |
471 | element = table->GetElement(13); | |
472 | mix->DefineElement(4, element, 0.68572); | |
473 | break; | |
474 | case 14: //SDD hybrid microcab | |
475 | mix = new TGeoMixture("HYB-CAB", 5, 2.0502); | |
476 | element = table->GetElement(1); | |
477 | mix->DefineElement(0, element, 0.00926228815); | |
478 | element = table->GetElement(6); | |
479 | mix->DefineElement(1, element, 0.24281879711); | |
480 | element = table->GetElement(7); | |
481 | mix->DefineElement(2, element, 0.02574224025); | |
482 | element = table->GetElement(8); | |
483 | mix->DefineElement(3, element, 0.07348667449); | |
484 | element = table->GetElement(13); | |
485 | mix->DefineElement(4, element, 0.64869); | |
486 | break; | |
487 | case 15: //SDD anode microcab | |
488 | mix = new TGeoMixture("ANOD-CAB", 5, 1.7854); | |
489 | element = table->GetElement(1); | |
490 | mix->DefineElement(0, element, 0.0128595919215); | |
491 | element = table->GetElement(6); | |
492 | mix->DefineElement(1, element, 0.392653705471); | |
493 | element = table->GetElement(7); | |
494 | mix->DefineElement(2, element, 0.041626868025); | |
495 | element = table->GetElement(8); | |
496 | mix->DefineElement(3, element, 0.118832707289); | |
497 | element = table->GetElement(13); | |
498 | mix->DefineElement(4, element, 0.431909); | |
499 | break; | |
500 | case 16: // inox/alum | |
501 | mix = new TGeoMixture("INOX-AL", 5, 3.0705); | |
502 | element = table->GetElement(13); | |
503 | mix->DefineElement(0, element, 0.816164); | |
504 | element = table->GetElement(14); | |
505 | mix->DefineElement(1, element, 0.000919182); | |
506 | element = table->GetElement(24); | |
507 | mix->DefineElement(2, element, 0.0330906); | |
508 | element = table->GetElement(26); | |
509 | mix->DefineElement(3, element, 0.131443); | |
510 | element = table->GetElement(28); | |
511 | mix->DefineElement(4, element, 0.0183836); | |
512 | case 17: // MYLAR | |
513 | mix = new TGeoMixture("MYLAR", 3, 1.39); | |
514 | element = table->GetElement(1); | |
515 | mix->DefineElement(0, element, 0.0416667); | |
516 | element = table->GetElement(6); | |
517 | mix->DefineElement(1, element, 0.625); | |
518 | element = table->GetElement(8); | |
519 | mix->DefineElement(2, element, 0.333333); | |
520 | break; | |
521 | case 18: // SPDBUS(AL+KPT+EPOX) - unknown composition | |
522 | mix = new TGeoMixture("SPDBUS", 1, 1.906); | |
523 | element = table->GetElement(9); | |
524 | mix->DefineElement(0, element, 1.); | |
525 | break; | |
526 | case 19: // SDD/SSD rings - unknown composition | |
527 | mix = new TGeoMixture("SDDRINGS", 1, 1.8097); | |
528 | element = table->GetElement(6); | |
529 | mix->DefineElement(0, element, 1.); | |
530 | break; | |
531 | case 20: // SPD end ladder - unknown composition | |
532 | mix = new TGeoMixture("SPDEL", 1, 3.6374); | |
533 | element = table->GetElement(22); | |
534 | mix->DefineElement(0, element, 1.); | |
535 | break; | |
536 | case 21: // SDD end ladder - unknown composition | |
537 | mix = new TGeoMixture("SDDEL", 1, 0.3824); | |
538 | element = table->GetElement(30); | |
539 | mix->DefineElement(0, element, 1.); | |
540 | break; | |
541 | case 22: // SSD end ladder - unknown composition | |
542 | mix = new TGeoMixture("SSDEL", 1, 0.68); | |
543 | element = table->GetElement(16); | |
544 | mix->DefineElement(0, element, 1.); | |
545 | break; | |
546 | } | |
547 | mix->SetZ(z); | |
548 | printf("Patched with mixture %s\n", mix->GetName()); | |
549 | return mix; | |
8495a208 | 550 | } |
551 | ||
8495a208 | 552 | //_____________________________________________________________________________ |
efde9b4d | 553 | void TFlukaMCGeometry::CreateFlukaMatFile(const char *fname) |
8495a208 | 554 | { |
555 | // ==== from FLUGG ==== | |
556 | // NAMES OF ELEMENTS AND COMPOUNDS: the names must be written in upper case, | |
557 | // according to the fluka standard. In addition,. they must be equal to the | |
558 | // names of the fluka materials - see fluka manual - in order that the | |
559 | // program load the right cross sections, and equal to the names included in | |
560 | // the .pemf. Otherwise the user must define the LOW-MAT CARDS, and make his | |
561 | // own .pemf, in order to get the right cross sections loaded in memory. | |
89125dc7 | 562 | |
563 | ||
8495a208 | 564 | TString sname; |
565 | gGeoManager->Export("flgeom.root"); | |
566 | if (fname) sname = fname; | |
567 | else sname = "flukaMat.inp"; | |
568 | ofstream out; | |
569 | out.open(sname.Data(), ios::out); | |
570 | if (!out.good()) { | |
571 | Fatal("CreateFlukaMatFile", "could not open file %s for writing", sname.Data()); | |
572 | return; | |
573 | } | |
574 | PrintHeader(out, "MATERIALS AND COMPOUNDS"); | |
575 | PrintHeader(out, "MATERIALS"); | |
89125dc7 | 576 | Int_t i,j,idmat; |
577 | Int_t counttothree, nelem; | |
578 | Double_t a,z,rho, w; | |
579 | TGeoElementTable *table = TGeoElementTable::Instance(); | |
580 | TGeoElement *element; | |
581 | element = table->GetElement(13); | |
582 | element->SetTitle("ALUMINUM"); // this is how FLUKA likes it ... | |
583 | element = table->GetElement(15); | |
584 | element->SetTitle("PHOSPHO"); // same story ... | |
585 | // element = table->GetElement(10); | |
586 | // element->SetTitle("ARGON"); // NEON not in neutron xsec table | |
587 | Int_t nelements = table->GetNelements(); | |
8495a208 | 588 | TList *matlist = gGeoManager->GetListOfMaterials(); |
89125dc7 | 589 | TList *medlist = gGeoManager->GetListOfMedia(); |
590 | Int_t nmed = medlist->GetSize(); | |
8495a208 | 591 | TIter next(matlist); |
592 | Int_t nmater = matlist->GetSize(); | |
593 | Int_t nfmater = 0; | |
89125dc7 | 594 | TGeoMaterial *mat; |
8495a208 | 595 | TGeoMixture *mix = 0; |
596 | TString matname; | |
89125dc7 | 597 | TObjString *objstr; |
598 | // Create all needed elements | |
599 | for (Int_t i=1; i<nelements; i++) { | |
600 | element = table->GetElement(i); | |
601 | // skip elements which are not defined | |
602 | if (!element->IsUsed() && !element->IsDefined()) continue; | |
603 | matname = element->GetTitle(); | |
604 | ToFlukaString(matname); | |
605 | rho = 0.999; | |
606 | mat = new TGeoMaterial(matname, element->A(), element->Z(), rho); | |
607 | mat->SetIndex(nfmater+3); | |
608 | mat->SetUsed(kTRUE); | |
609 | fMatList->Add(mat); | |
610 | objstr = new TObjString(matname.Data()); | |
611 | fMatNames->Add(objstr); | |
612 | nfmater++; | |
613 | } | |
614 | Int_t indmat = nfmater; | |
615 | TGeoMedium *med; | |
616 | // Adjust material names and add them to FLUKA list | |
8495a208 | 617 | for (i=0; i<nmater; i++) { |
618 | mat = (TGeoMaterial*)matlist->At(i); | |
89125dc7 | 619 | if (!mat->IsUsed()) continue; |
620 | z = mat->GetZ(); | |
621 | a = mat->GetA(); | |
622 | rho = mat->GetDensity(); | |
623 | if (mat->GetZ()<0.001) { | |
efde9b4d | 624 | mat->SetIndex(2); // vacuum, built-in inside FLUKA |
625 | continue; | |
89125dc7 | 626 | } |
627 | matname = mat->GetName(); | |
628 | FlukaMatName(matname); | |
629 | // material with one element: create it as mixture since it can be duplicated | |
630 | if (!mat->IsMixture()) { | |
631 | // normal materials | |
632 | mix = new TGeoMixture(matname.Data(), 1, rho); | |
633 | mix->DefineElement(0, mat->GetElement(), 1.); | |
8495a208 | 634 | mat->SetIndex(nfmater+3); |
89125dc7 | 635 | for (j=0; j<nmed; j++) { |
636 | med = (TGeoMedium*)medlist->At(j); | |
637 | if (med->GetMaterial() == mat) { | |
638 | med->SetMaterial(mix); | |
639 | if (mat->GetCerenkovProperties()) { | |
640 | mix->SetCerenkovProperties(mat->GetCerenkovProperties()); | |
641 | mat->SetCerenkovProperties(0); | |
8495a208 | 642 | } |
89125dc7 | 643 | break; |
644 | } | |
645 | } | |
646 | mat = (TGeoMaterial*)mix; | |
8495a208 | 647 | } |
89125dc7 | 648 | mat->SetIndex(nfmater+3); |
649 | objstr = new TObjString(matname.Data()); | |
650 | fMatList->Add(mat); | |
651 | fMatNames->Add(objstr); | |
652 | nfmater++; | |
653 | } | |
654 | ||
655 | // Dump all elements with MATERIAL cards | |
8495a208 | 656 | for (i=0; i<nfmater; i++) { |
89125dc7 | 657 | mat = (TGeoMaterial*)fMatList->At(i); |
658 | // mat->SetUsed(kFALSE); | |
659 | mix = 0; | |
660 | // out << "* " << mat->GetName() << endl; | |
8495a208 | 661 | out << setw(10) << "MATERIAL "; |
662 | out.setf(static_cast<std::ios::fmtflags>(0),std::ios::floatfield); | |
89125dc7 | 663 | objstr = (TObjString*)fMatNames->At(i); |
b0853588 | 664 | matname = objstr->GetString(); |
89125dc7 | 665 | z = mat->GetZ(); |
666 | a = mat->GetA(); | |
667 | rho = mat->GetDensity(); | |
8495a208 | 668 | if (mat->IsMixture()) { |
669 | out << setw(10) << " "; | |
670 | out << setw(10) << " "; | |
671 | mix = (TGeoMixture*)mat; | |
672 | } else { | |
89125dc7 | 673 | out << setw(10) << setiosflags(ios::fixed) << setprecision(1) << z; |
674 | out << setw(10) << setprecision(3) << a; | |
8495a208 | 675 | } |
676 | out.setf(static_cast<std::ios::fmtflags>(0),std::ios::floatfield); | |
89125dc7 | 677 | out << setw(10) << setiosflags(ios::scientific) << setprecision(3) << rho; |
8495a208 | 678 | out.setf(static_cast<std::ios::fmtflags>(0),std::ios::floatfield); |
89125dc7 | 679 | out << setw(10) << setiosflags(ios::fixed) << setprecision(1) << Double_t(mat->GetIndex()); |
8495a208 | 680 | out << setw(10) << " "; |
efde9b4d | 681 | out << setw(10) << " "; |
682 | out << setw(8) << matname.Data() << endl; | |
89125dc7 | 683 | if (!mix) { |
684 | // add LOW-MAT card for NEON to associate with ARGON neutron xsec | |
685 | if (z==10 && matname.Contains("NEON")) { | |
686 | out << setw(10) << "LOW-MAT "; | |
687 | out.setf(static_cast<std::ios::fmtflags>(0),std::ios::floatfield); | |
688 | out << setw(10) << setiosflags(ios::fixed) << setprecision(1) << Double_t(mat->GetIndex()); | |
689 | out << setw(10) << setiosflags(ios::fixed) << setprecision(1) << 18.; | |
690 | out << setw(10) << setiosflags(ios::fixed) << setprecision(1) << -2.; | |
691 | out << setw(10) << setiosflags(ios::fixed) << setprecision(1) << 293.; | |
692 | out << setw(10) << " "; | |
693 | out << setw(10) << " "; | |
694 | // out << setw(8) << matname.Data() << endl; | |
695 | out << setw(8) << " " << endl; | |
696 | } | |
697 | continue; | |
d11bf3ca | 698 | } |
8495a208 | 699 | counttothree = 0; |
700 | out << setw(10) << "COMPOUND "; | |
701 | nelem = mix->GetNelements(); | |
89125dc7 | 702 | objstr = (TObjString*)fMatNames->At(i); |
8495a208 | 703 | matname = objstr->GetString(); |
8495a208 | 704 | for (j=0; j<nelem; j++) { |
89125dc7 | 705 | w = (mix->GetWmixt())[j]; |
706 | if (w<0.00001) w=0.00001; | |
707 | z = (mix->GetZmixt())[j]; | |
708 | a = (mix->GetAmixt())[j]; | |
709 | idmat = GetElementIndex(Int_t(z)); | |
710 | if (!idmat) Error("CreateFlukaMatFile", "element with Z=%f not found", z); | |
8495a208 | 711 | out.setf(static_cast<std::ios::fmtflags>(0),std::ios::floatfield); |
89125dc7 | 712 | out << setw(10) << setiosflags(ios::fixed) << setprecision(6) << -w; |
8495a208 | 713 | out.setf(static_cast<std::ios::fmtflags>(0),std::ios::floatfield); |
714 | out << setw(10) << setiosflags(ios::fixed) << setprecision(1) << Double_t(idmat); | |
715 | counttothree++; | |
716 | if (counttothree == 3) { | |
717 | out << matname.Data(); | |
718 | out << endl; | |
719 | if ( (j+1) != nelem) out << setw(10) << "COMPOUND "; | |
720 | counttothree = 0; | |
89125dc7 | 721 | } |
722 | } | |
8495a208 | 723 | if (nelem%3) { |
724 | for (j=0; j<(3-(nelem%3)); j++) | |
725 | out << setw(10) << " " << setw(10) << " "; | |
726 | out << matname.Data(); | |
727 | out << endl; | |
89125dc7 | 728 | } |
729 | } | |
8495a208 | 730 | Int_t nvols = gGeoManager->GetListOfUVolumes()->GetEntriesFast()-1; |
731 | TGeoVolume *vol; | |
8495a208 | 732 | // Now print the material assignments |
d11bf3ca | 733 | Double_t flagfield = 0.; |
734 | printf("#############################################################\n"); | |
89125dc7 | 735 | if (gFluka->IsFieldEnabled()) { |
d11bf3ca | 736 | flagfield = 1.; |
737 | printf("Magnetic field enabled\n"); | |
738 | } else printf("Magnetic field disabled\n"); | |
739 | printf("#############################################################\n"); | |
740 | ||
8495a208 | 741 | PrintHeader(out, "TGEO MATERIAL ASSIGNMENTS"); |
742 | for (i=1; i<=nvols; i++) { | |
743 | vol = gGeoManager->GetVolume(i); | |
744 | mat = vol->GetMedium()->GetMaterial(); | |
89125dc7 | 745 | // mat->SetUsed(kTRUE); |
8495a208 | 746 | idmat = mat->GetIndex(); |
89125dc7 | 747 | for (Int_t j=0; j<nfmater; j++) { |
748 | mat = (TGeoMaterial*)fMatList->At(j); | |
749 | if (mat->GetIndex() == idmat) mat->SetUsed(kTRUE); | |
750 | } | |
8495a208 | 751 | out << setw(10) << "ASSIGNMAT "; |
752 | out.setf(static_cast<std::ios::fmtflags>(0),std::ios::floatfield); | |
753 | out << setw(10) << setiosflags(ios::fixed) << Double_t(idmat); | |
754 | out << setw(10) << setiosflags(ios::fixed) << Double_t(i); | |
755 | out << setw(10) << "0.0"; | |
b0853588 | 756 | out << setw(10) << "0.0"; |
8495a208 | 757 | out << setw(10) << setiosflags(ios::fixed) << flagfield; |
b0853588 | 758 | out << setw(10) << "0.0"; |
8495a208 | 759 | out << endl; |
760 | } | |
8495a208 | 761 | out.close(); |
efde9b4d | 762 | fLastMaterial = nfmater+2; |
89125dc7 | 763 | /* |
764 | TGeoMaterial *mat1 = 0; | |
765 | for (i=1; i<=nvols; i++) { | |
766 | vol = gGeoManager->GetVolume(i); | |
767 | med = vol->GetMedium(); | |
768 | mat = med->GetMaterial(); | |
769 | printf("Region %d: %s\n", i, vol->GetName()); | |
770 | printf(" medium %d: %s\n", med->GetId(), med->GetName()); | |
771 | for (j=0; j<nfmater; j++) { | |
772 | mat1 = (TGeoMaterial*)fMatList->At(j); | |
773 | if (mat1 != mat) continue; | |
774 | objstr = (TObjString*)fMatNames->At(j); | |
775 | matname = objstr->GetString(); | |
776 | break; | |
777 | } | |
778 | if (mat1 != mat) printf(" (*) material not found in Fluka list\n"); | |
779 | printf(" material %s (at ind=%d): FlukaID=%d FlukaName=%s\n", | |
780 | mat->GetName(), j, mat->GetIndex(), matname.Data()); | |
781 | if (mat->GetCerenkovProperties()) printf(" Cerenkov properties found\n"); | |
782 | } | |
783 | */ | |
784 | ||
785 | if (!gFluka->IsGeneratePemf()) return; | |
786 | // Write peg files | |
787 | char number[20]; | |
788 | for (i=indmat; i<nfmater; i++) { | |
789 | mat = (TGeoMaterial*)fMatList->At(i); | |
790 | if (!mat->IsUsed()) continue; | |
791 | sname = "mat"; | |
792 | sprintf(number, "%d", i); | |
793 | sname.Append(number); | |
794 | WritePegFile(i); | |
795 | sname.Prepend("$FLUPRO/pemf/rpemf peg/"); | |
796 | gSystem->Exec(sname.Data()); | |
797 | } | |
798 | sname = "cat peg/*.pemf > peg/alice.pemf"; | |
799 | gSystem->Exec(sname.Data()); | |
800 | sname = "mv peg/alice.pemf alice.pemf"; | |
801 | gSystem->Exec(sname.Data()); | |
802 | } | |
803 | ||
804 | //_____________________________________________________________________________ | |
805 | void TFlukaMCGeometry::WritePegFile(Int_t imat) const | |
806 | { | |
807 | // Write the .peg file for one material | |
808 | TGeoMaterial *mat = (TGeoMaterial*)fMatList->At(imat); | |
809 | TString name = ((TObjString*)fMatNames->At(imat))->GetString(); | |
810 | TString line; | |
811 | char number[20]; | |
812 | TGeoElement *elem = mat->GetElement(); | |
813 | name = name.Strip(); | |
814 | TString sname = "mat"; | |
815 | sprintf(number, "%d", imat); | |
816 | sname.Append(number); | |
817 | sname.Append(".peg"); | |
818 | sname.Prepend("peg/"); | |
819 | ofstream out; | |
820 | out.open(sname.Data(), ios::out); | |
821 | if (!out.good()) return; | |
822 | Double_t dens = mat->GetDensity(); | |
823 | TGeoMixture *mix = 0; | |
824 | Int_t nel = 1; | |
825 | Int_t i; | |
826 | if (mat->IsMixture()) { | |
827 | mix = (TGeoMixture*)mat; | |
828 | nel = mix->GetNelements(); | |
829 | } | |
830 | if (nel==1) { | |
831 | out << "ELEM" << endl; | |
832 | out << " &INP IRAYL=1, RHO=" << dens << ", " << endl; | |
833 | if (dens<0.01) out << " GASP=1." << endl; | |
834 | out << " &END" << endl; | |
835 | out << name.Data() << endl; | |
836 | out << elem->GetName() << endl; | |
837 | } else { | |
838 | out << "MIXT" << endl; | |
839 | out << " &INP IRAYL=1, NE=" << nel << ", RHOZ="; | |
840 | line = ""; | |
841 | for (i=0; i<nel; i++) { | |
842 | sprintf(number, "%f", mix->GetWmixt()[i]); | |
843 | line += number; | |
844 | line += ", "; | |
845 | if (line.Length() > 30) { | |
846 | out << line.Data() << endl; | |
847 | line = ""; | |
848 | } | |
849 | } | |
850 | if (line.Length()) out << " " << line.Data() << endl; | |
851 | out << " RHO=" << dens; | |
852 | if (dens<0.01) out << ", GASP=1." << endl; | |
853 | out << " &END" << endl; | |
854 | out << name.Data() << endl; | |
855 | for (i=0; i<nel; i++) { | |
856 | elem = mix->GetElement(i); | |
857 | line = elem->GetName(); | |
858 | if (line.Length()==1) line.Append(" "); | |
859 | out << line.Data() << " "; | |
860 | } | |
861 | out << endl; | |
862 | } | |
863 | out << "ENER" << endl; | |
864 | out << " $INP AE=0.56099906, UE=3000000., AP=.03, UP=3000000. $END" << endl; | |
865 | out << "PWLF" << endl; | |
866 | out << " $INP NALE=300, NALG=400, NALR=100 $END" << endl; | |
867 | out << "DECK" << endl; | |
868 | out << " $INP $END" << endl; | |
869 | out << "TEST" << endl; | |
870 | out << " $INP $END" << endl; | |
871 | out.close(); | |
8495a208 | 872 | } |
873 | ||
874 | //_____________________________________________________________________________ | |
875 | void TFlukaMCGeometry::PrintHeader(ofstream &out, const char *text) const | |
876 | { | |
877 | // Print a FLUKA header. | |
878 | out << "*\n" << "*\n" << "*\n"; | |
879 | out << "********************* " << text << " *********************\n" | |
880 | << "*\n"; | |
881 | out << "*...+....1....+....2....+....3....+....4....+....5....+....6....+....7..." | |
882 | << endl; | |
883 | out << "*" << endl; | |
884 | } | |
885 | ||
886 | //_____________________________________________________________________________ | |
887 | Int_t TFlukaMCGeometry::RegionId() const | |
888 | { | |
889 | // Returns current region id <-> TGeo node id | |
890 | if (gGeoManager->IsOutside()) return 0; | |
891 | return gGeoManager->GetCurrentNode()->GetUniqueID(); | |
892 | } | |
89125dc7 | 893 | |
894 | //_____________________________________________________________________________ | |
895 | Int_t TFlukaMCGeometry::GetElementIndex(Int_t z) const | |
896 | { | |
6f5667d1 | 897 | // Get index of a material having a given Z element. |
89125dc7 | 898 | TIter next(fMatList); |
899 | TGeoMaterial *mat; | |
900 | Int_t index = 0; | |
901 | while ((mat=(TGeoMaterial*)next())) { | |
902 | if (mat->IsMixture()) continue; | |
903 | if (mat->GetElement()->Z() == z) return mat->GetIndex(); | |
904 | } | |
905 | return index; | |
906 | } | |
907 | ||
05265ca9 | 908 | //_____________________________________________________________________________ |
909 | void TFlukaMCGeometry::SetMreg(Int_t mreg) | |
910 | { | |
911 | // Update if needed next history; | |
89125dc7 | 912 | if (gFluka->GetDummyBoundary()==2) { |
d11bf3ca | 913 | gGeoManager->CdNode(fNextLattice-1); |
914 | return; | |
915 | } | |
89125dc7 | 916 | Int_t curreg = (gGeoManager->IsOutside())?(gMCGeom->NofVolumes()+1):gGeoManager->GetCurrentVolume()->GetNumber(); |
05265ca9 | 917 | if (mreg==curreg) return; |
918 | if (mreg==fNextRegion) { | |
919 | if (fNextLattice!=999999999) gGeoManager->CdNode(fNextLattice-1); | |
920 | return; | |
d11bf3ca | 921 | } else { |
922 | if (mreg == fCurrentRegion) { | |
923 | if (fCurrentLattice!=999999999) gGeoManager->CdNode(fCurrentLattice-1); | |
924 | return; | |
925 | } | |
926 | } | |
2bc4c610 | 927 | if (fDebug) printf("ERROR: mreg=%i neither current nor next region\n", mreg); |
05265ca9 | 928 | } |
929 | ||
d11bf3ca | 930 | //_____________________________________________________________________________ |
931 | void TFlukaMCGeometry::SetCurrentRegion(Int_t mreg, Int_t latt) | |
932 | { | |
933 | // Set index/history for next entered region | |
934 | fCurrentRegion = mreg; | |
935 | fCurrentLattice = latt; | |
936 | } | |
937 | ||
05265ca9 | 938 | //_____________________________________________________________________________ |
939 | void TFlukaMCGeometry::SetNextRegion(Int_t mreg, Int_t latt) | |
940 | { | |
941 | // Set index/history for next entered region | |
942 | fNextRegion = mreg; | |
943 | fNextLattice = latt; | |
944 | } | |
8495a208 | 945 | |
946 | //_____________________________________________________________________________ | |
947 | void TFlukaMCGeometry::ToFlukaString(TString &str) const | |
948 | { | |
949 | // ToFlukaString converts an string to something usefull in FLUKA: | |
950 | // * Capital letters | |
951 | // * Only 8 letters | |
952 | // * Replace ' ' by '_' | |
953 | if (str.Length()<8) { | |
954 | str += " "; | |
955 | } | |
956 | str.Remove(8); | |
957 | Int_t ilast; | |
958 | for (ilast=7; ilast>0; ilast--) if (str(ilast)!=' ') break; | |
959 | str.ToUpper(); | |
960 | for (Int_t pos=0; pos<ilast; pos++) | |
961 | if (str(pos)==' ') str.Replace(pos,1,"_",1); | |
962 | return; | |
963 | } | |
89125dc7 | 964 | |
965 | //_____________________________________________________________________________ | |
966 | void TFlukaMCGeometry::FlukaMatName(TString &str) const | |
967 | { | |
6f5667d1 | 968 | // Convert a name to upper case 8 chars. |
89125dc7 | 969 | ToFlukaString(str); |
970 | Int_t ilast; | |
971 | for (ilast=7; ilast>0; ilast--) if (str(ilast)!=' ') break; | |
972 | if (ilast>5) ilast = 5; | |
973 | char number[3]; | |
974 | TIter next(fMatNames); | |
975 | TObjString *objstr; | |
976 | TString matname; | |
977 | Int_t index = 0; | |
978 | while ((objstr=(TObjString*)next())) { | |
979 | matname = objstr->GetString(); | |
980 | if (matname == str) { | |
981 | index++; | |
982 | if (index<10) { | |
983 | number[0] = '0'; | |
984 | sprintf(&number[1], "%d", index); | |
985 | } else if (index<100) { | |
986 | sprintf(number, "%d", index); | |
987 | } else { | |
988 | Error("FlukaMatName", "Too many materials %s", str.Data()); | |
989 | return; | |
990 | } | |
991 | str.Replace(ilast+1, 2, number); | |
992 | str.Remove(8); | |
993 | } | |
994 | } | |
995 | } | |
996 | ||
8495a208 | 997 | //______________________________________________________________________________ |
998 | void TFlukaMCGeometry::Vname(const char *name, char *vname) const | |
999 | { | |
1000 | // | |
1001 | // convert name to upper case. Make vname at least 4 chars | |
1002 | // | |
1003 | Int_t l = strlen(name); | |
1004 | Int_t i; | |
1005 | l = l < 4 ? l : 4; | |
1006 | for (i=0;i<l;i++) vname[i] = toupper(name[i]); | |
1007 | for (i=l;i<4;i++) vname[i] = ' '; | |
1008 | vname[4] = 0; | |
1009 | } | |
1010 | ||
1011 | ||
1012 | // FLUKA GEOMETRY WRAPPERS - to replace FLUGG wrappers | |
1013 | ||
8495a208 | 1014 | //_____________________________________________________________________________ |
1015 | Int_t idnrwr(const Int_t & /*nreg*/, const Int_t & /*mlat*/) | |
1016 | { | |
1017 | // from FLUGG: | |
1018 | // Wrapper for setting DNEAR option on fluka side. Must return 0 | |
1019 | // if user doesn't want Fluka to use DNEAR to compute the | |
1020 | // step (the same effect is obtained with the GLOBAL (WHAT(3)=-1) | |
1021 | // card in fluka input), returns 1 if user wants Fluka always to | |
1022 | // use DNEAR (in this case, be sure that GEANT4 DNEAR is unique, | |
1023 | // coming from all directions!!!) | |
89125dc7 | 1024 | if (gMCGeom->IsDebugging()) printf("========== Dummy IDNRWR\n"); |
8495a208 | 1025 | return 0; |
1026 | } | |
1027 | ||
8495a208 | 1028 | //_____________________________________________________________________________ |
1029 | void g1wr(Double_t &pSx, Double_t &pSy, Double_t &pSz, | |
d11bf3ca | 1030 | Double_t *pV, Int_t &oldReg , const Int_t &oldLttc, Double_t &propStep, |
1031 | Int_t &nascFlag, Double_t &retStep, Int_t &newReg, | |
2573ac89 | 1032 | Double_t &saf, Int_t &newLttc, Int_t <tcFlag, |
8495a208 | 1033 | Double_t *sLt, Int_t *jrLt) |
d11bf3ca | 1034 | |
8495a208 | 1035 | { |
d11bf3ca | 1036 | // Initialize FLUKa point and direction; |
6f5667d1 | 1037 | gNstep++; |
d11bf3ca | 1038 | /* |
1039 | if (kNstep>0) { | |
89125dc7 | 1040 | gMCGeom->SetDebugMode(kTRUE); |
1041 | gFluka->SetVerbosityLevel(3); | |
d11bf3ca | 1042 | } |
1043 | if (kNstep>6520) { | |
89125dc7 | 1044 | gMCGeom->SetDebugMode(kFALSE); |
1045 | gFluka->SetVerbosityLevel(0); | |
d11bf3ca | 1046 | } |
1047 | if ((kNstep%10)==0) printf("step %i\n", kNstep); | |
1048 | */ | |
8495a208 | 1049 | |
89125dc7 | 1050 | if (gMCGeom->IsDebugging()) { |
2bc4c610 | 1051 | printf("========== Inside G1WR\n"); |
1052 | printf(" point/dir:(%14.9f, %14.9f, %14.9f, %g, %g, %g)\n", pSx,pSy,pSz,pV[0],pV[1],pV[2]); | |
d11bf3ca | 1053 | printf(" oldReg=%i oldLttc=%i pstep=%f\n",oldReg, oldLttc, propStep); |
2bc4c610 | 1054 | } |
8495a208 | 1055 | gGeoManager->SetCurrentPoint(pSx, pSy, pSz); |
1056 | gGeoManager->SetCurrentDirection(pV); | |
89125dc7 | 1057 | gMCGeom->SetCurrentRegion(oldReg, oldLttc); |
d11bf3ca | 1058 | // Initialize default return values |
1059 | lttcFlag = 0; | |
1060 | jrLt[lttcFlag] = oldLttc; | |
1061 | sLt[lttcFlag] = propStep; | |
1062 | jrLt[lttcFlag+1] = -1; | |
1063 | sLt[lttcFlag+1] = 0.; | |
1064 | newReg = oldReg; | |
1065 | newLttc = oldLttc; | |
1066 | // check if dummy boundary flag is set | |
1067 | Int_t curLttc, curReg; | |
89125dc7 | 1068 | if (gFluka->IsDummyBoundary()) { |
d11bf3ca | 1069 | // printf("Dummy boundary intercepted. Point is: %f, %f, %f\n", pSx, pSy, pSz); |
1070 | Bool_t crossedDummy = (oldLttc == TFlukaMCGeometry::kLttcVirtual)?kTRUE:kFALSE; | |
1071 | if (crossedDummy) { | |
1072 | // FLUKA crossed the dummy boundary - update new region/history | |
1073 | retStep = 0.; | |
1074 | saf = 0.; | |
89125dc7 | 1075 | gMCGeom->GetNextRegion(newReg, newLttc); |
1076 | gMCGeom->SetMreg(newReg); | |
1077 | if (gMCGeom->IsDebugging()) printf(" virtual newReg=%i newLttc=%i\n", newReg, newLttc); | |
d11bf3ca | 1078 | sLt[lttcFlag] = 0.; // null step in current region |
1079 | lttcFlag++; | |
1080 | jrLt[lttcFlag] = newLttc; | |
1081 | sLt[lttcFlag] = 0.; // null step in next region | |
1082 | jrLt[lttcFlag+1] = -1; | |
1083 | sLt[lttcFlag+1] = 0.; | |
89125dc7 | 1084 | gFluka->SetDummyBoundary(0); |
d11bf3ca | 1085 | return; |
1086 | } | |
1087 | } | |
1088 | ||
1089 | // Reset outside flag | |
05265ca9 | 1090 | if (gGeoManager->IsOutside()) { |
1091 | gGeoManager->SetOutside(kFALSE); | |
1092 | gGeoManager->CdTop(); | |
d11bf3ca | 1093 | } |
1094 | ||
1095 | // Reset dummy boundary flag | |
89125dc7 | 1096 | gFluka->SetDummyBoundary(0); |
d11bf3ca | 1097 | |
1098 | curLttc = gGeoManager->GetCurrentNodeId()+1; | |
1099 | curReg = gGeoManager->GetCurrentVolume()->GetNumber(); | |
2573ac89 | 1100 | if (oldLttc != curLttc) { |
d11bf3ca | 1101 | // FLUKA crossed the boundary : we trust that the given point is really there, |
1102 | // so we just update TGeo state | |
2573ac89 | 1103 | gGeoManager->CdNode(oldLttc-1); |
1104 | curLttc = gGeoManager->GetCurrentNodeId()+1; | |
d11bf3ca | 1105 | curReg = gGeoManager->GetCurrentVolume()->GetNumber(); |
89125dc7 | 1106 | if (gMCGeom->IsDebugging()) printf(" re-initialized point: curReg=%i curLttc=%i\n", curReg, curLttc); |
d11bf3ca | 1107 | } |
1108 | // Now the current TGeo state reflects the FLUKA state | |
89125dc7 | 1109 | if (gMCGeom->IsDebugging()) printf(" current path: %s\n", gGeoManager->GetPath()); |
d11bf3ca | 1110 | Double_t extra = 1E-6; |
1111 | Double_t tmpStep = propStep + extra; | |
1112 | gGeoManager->FindNextBoundary(-tmpStep); | |
1113 | Double_t snext = gGeoManager->GetStep(); | |
1114 | // !!!!! | |
1115 | if (snext<=0) { | |
1116 | // FLUKA is in the wrong region, notify it | |
89125dc7 | 1117 | if (gMCGeom->IsDebugging()) printf("ERROR: snext=%f\n", snext); |
d11bf3ca | 1118 | // newReg = -3; |
1119 | // return; | |
1120 | snext = extra; | |
1121 | } | |
1122 | saf = gGeoManager->GetSafeDistance(); | |
1123 | Bool_t cross = kFALSE; | |
1124 | Bool_t onBound = kFALSE; | |
1125 | if (snext<tmpStep) { | |
1126 | // We have some boundary in the way | |
1127 | Double_t dd = snext-propStep; | |
1128 | if (dd < 0) { | |
1129 | cross = kTRUE; | |
1130 | dd = -dd; | |
1131 | } | |
1132 | if (dd < 1E-8) onBound = kTRUE; | |
1133 | } | |
1134 | snext += 1.E-8; | |
89125dc7 | 1135 | if (gMCGeom->IsDebugging()) { |
d11bf3ca | 1136 | if (!cross) printf(" physical step approved: %f\n", propStep); |
1137 | else printf(" boundary crossing at: %f\n", snext); | |
1138 | if (onBound) printf(" step on boundary limit ! NASC=%i\n", nascFlag); | |
1139 | } | |
1140 | if (!cross) { | |
1141 | // Next boundary further than proposed step, which is approved | |
1142 | retStep = propStep; | |
1143 | sLt[lttcFlag] = propStep; | |
1144 | return; | |
1145 | } | |
1146 | // The next boundary is closer. We try to cross it. | |
2573ac89 | 1147 | Double_t *point = gGeoManager->GetCurrentPoint(); |
1148 | Double_t *dir = gGeoManager->GetCurrentDirection(); | |
d11bf3ca | 1149 | Double_t pt[3]; |
1150 | memcpy(pt, point, 3*sizeof(Double_t)); | |
1151 | ||
2573ac89 | 1152 | Int_t i; |
d11bf3ca | 1153 | for (i=0;i<3;i++) point[i] += snext*dir[i]; |
1154 | gGeoManager->FindNode(); | |
1155 | newLttc = (gGeoManager->IsOutside())?(TFlukaMCGeometry::kLttcOutside):gGeoManager->GetCurrentNodeId()+1; | |
1156 | if (newLttc == oldLttc) { | |
1157 | // brute force ... | |
1158 | // Just try a fast extra small step | |
1159 | snext += 1E-6; | |
1160 | for (i=0;i<3;i++) point[i] = pt[i]+snext*dir[i]; | |
1161 | gGeoManager->FindNode(); | |
1162 | newLttc = (gGeoManager->IsOutside())?(TFlukaMCGeometry::kLttcOutside):gGeoManager->GetCurrentNodeId()+1; | |
1163 | if (newLttc == oldLttc) { | |
1164 | // check if state changes at the end of the proposed step | |
1165 | for (i=0;i<3;i++) point[i] = pt[i]+propStep*dir[i]; | |
2573ac89 | 1166 | gGeoManager->FindNode(); |
d11bf3ca | 1167 | newLttc = (gGeoManager->IsOutside())?(TFlukaMCGeometry::kLttcOutside):gGeoManager->GetCurrentNodeId()+1; |
1168 | if (newLttc==oldLttc) { | |
1169 | // approve step | |
1170 | retStep = propStep; | |
1171 | sLt[lttcFlag] = propStep; | |
1172 | return; | |
05265ca9 | 1173 | } |
d11bf3ca | 1174 | // snext is underestimated - we will create a virtual one to overcome the error |
1175 | // printf("some boundary in the way...\n"); | |
1176 | } | |
1177 | } | |
1178 | gGeoManager->SetCurrentPoint(pt); | |
1179 | // newLttc = (gGeoManager->IsOutside())?(TFlukaMCGeometry::kLttcOutside):gGeoManager->GetCurrentNodeId()+1; | |
89125dc7 | 1180 | newReg = (gGeoManager->IsOutside())?(gMCGeom->NofVolumes()+1):gGeoManager->GetCurrentVolume()->GetNumber(); |
1181 | if (gMCGeom->IsDebugging()) printf(" newReg=%i newLttc=%i\n", newReg, newLttc); | |
d11bf3ca | 1182 | |
1183 | // We really crossed the boundary, but is it the same region ? | |
89125dc7 | 1184 | gMCGeom->SetNextRegion(newReg, newLttc); |
d11bf3ca | 1185 | if (newReg == oldReg) { |
1186 | // Virtual boundary between replicants | |
89125dc7 | 1187 | if (gMCGeom->IsDebugging()) printf(" DUMMY boundary\n"); |
d11bf3ca | 1188 | newReg = 1; // cheat FLUKA telling it it crossed the TOP region |
1189 | newLttc = TFlukaMCGeometry::kLttcVirtual; | |
1190 | // mark that next boundary is virtual | |
89125dc7 | 1191 | gFluka->SetDummyBoundary(1); |
d11bf3ca | 1192 | } |
1193 | retStep = snext; | |
1194 | sLt[lttcFlag] = snext; | |
1195 | lttcFlag++; | |
1196 | jrLt[lttcFlag] = newLttc; | |
1197 | sLt[lttcFlag] = snext; | |
1198 | jrLt[lttcFlag+1] = -1; | |
1199 | sLt[lttcFlag+1] = 0.; | |
1200 | ||
2573ac89 | 1201 | if (newLttc!=oldLttc) { |
1202 | if (gGeoManager->IsOutside()) { | |
1203 | gGeoManager->SetOutside(kFALSE); | |
1204 | gGeoManager->CdTop(); | |
d11bf3ca | 1205 | } |
1206 | gGeoManager->CdTop(); | |
6f5667d1 | 1207 | if (!gGeoManager->GetCurrentMatrix()->IsIdentity()) printf("ERROR at step %i\n", gNstep); |
2573ac89 | 1208 | gGeoManager->CdNode(oldLttc-1); |
8495a208 | 1209 | } |
89125dc7 | 1210 | if (gMCGeom->IsDebugging()) { |
d11bf3ca | 1211 | printf("=> snext=%g safe=%g\n", snext, saf); |
1212 | for (Int_t i=0; i<lttcFlag+1; i++) printf(" jrLt[%i]=%i sLt[%i]=%g\n", i,jrLt[i],i,sLt[i]); | |
1213 | } | |
89125dc7 | 1214 | if (gMCGeom->IsDebugging()) printf("<= G1WR (in: %s)\n", gGeoManager->GetPath()); |
8495a208 | 1215 | } |
1216 | ||
efde9b4d | 1217 | //_____________________________________________________________________________ |
1218 | void g1rtwr() | |
1219 | { | |
89125dc7 | 1220 | if (gMCGeom->IsDebugging()) printf("========== Dummy G1RTWR\n"); |
efde9b4d | 1221 | } |
1222 | ||
1223 | //_____________________________________________________________________________ | |
1224 | void conhwr(Int_t & /*intHist*/, Int_t * /*incrCount*/) | |
1225 | { | |
89125dc7 | 1226 | if (gMCGeom->IsDebugging()) printf("========== Dummy CONHWR\n"); |
efde9b4d | 1227 | } |
1228 | ||
1229 | //_____________________________________________________________________________ | |
2573ac89 | 1230 | void inihwr(Int_t &intHist) |
efde9b4d | 1231 | { |
89125dc7 | 1232 | if (gMCGeom->IsDebugging()) printf("========== Inside INIHWR -> reinitializing history: %i\n", intHist); |
2573ac89 | 1233 | if (gGeoManager->IsOutside()) gGeoManager->CdTop(); |
1234 | if (intHist<=0) { | |
1235 | // printf("=== wrong history number\n"); | |
1236 | return; | |
1237 | } | |
1238 | if (intHist==0) gGeoManager->CdTop(); | |
1239 | else gGeoManager->CdNode(intHist-1); | |
89125dc7 | 1240 | if (gMCGeom->IsDebugging()) { |
2bc4c610 | 1241 | printf(" --- current path: %s\n", gGeoManager->GetPath()); |
1242 | printf("<= INIHWR\n"); | |
1243 | } | |
efde9b4d | 1244 | } |
1245 | ||
1246 | //_____________________________________________________________________________ | |
1247 | void jomiwr(const Int_t & /*nge*/, const Int_t & /*lin*/, const Int_t & /*lou*/, | |
1248 | Int_t &flukaReg) | |
1249 | { | |
1250 | // Geometry initialization wrapper called by FLUKAM. Provides to FLUKA the | |
1251 | // number of regions (volumes in TGeo) | |
1252 | // build application geometry | |
89125dc7 | 1253 | if (gMCGeom->IsDebugging()) printf("========== Inside JOMIWR\n"); |
2573ac89 | 1254 | flukaReg = gGeoManager->GetListOfUVolumes()->GetEntriesFast(); |
89125dc7 | 1255 | if (gMCGeom->IsDebugging()) printf("<= JOMIWR: last region=%i\n", flukaReg); |
efde9b4d | 1256 | } |
1257 | ||
1258 | //_____________________________________________________________________________ | |
1259 | void lkdbwr(Double_t &pSx, Double_t &pSy, Double_t &pSz, | |
1260 | Double_t * /*pV*/, const Int_t &oldReg, const Int_t &oldLttc, | |
1261 | Int_t &newReg, Int_t &flagErr, Int_t &newLttc) | |
1262 | { | |
89125dc7 | 1263 | if (gMCGeom->IsDebugging()) { |
2bc4c610 | 1264 | printf("========== Inside LKDBWR (%f, %f, %f)\n",pSx, pSy, pSz); |
1265 | // printf(" in: pV=(%f, %f, %f)\n", pV[0], pV[1], pV[2]); | |
1266 | printf(" in: oldReg=%i oldLttc=%i\n", oldReg, oldLttc); | |
1267 | } | |
efde9b4d | 1268 | TGeoNode *node = gGeoManager->FindNode(pSx, pSy, pSz); |
1269 | if (gGeoManager->IsOutside()) { | |
89125dc7 | 1270 | newReg = gMCGeom->NofVolumes()+1; |
05265ca9 | 1271 | // newLttc = gGeoManager->GetCurrentNodeId(); |
1272 | newLttc = 999999999; | |
89125dc7 | 1273 | if (gMCGeom->IsDebugging()) { |
2bc4c610 | 1274 | printf("OUTSIDE\n"); |
1275 | printf(" out: newReg=%i newLttc=%i\n", newReg, newLttc); | |
1276 | printf("<= LKMGWR\n"); | |
1277 | } | |
05265ca9 | 1278 | flagErr = newReg; |
1279 | return; | |
efde9b4d | 1280 | } |
1281 | newReg = node->GetVolume()->GetNumber(); | |
1282 | newLttc = gGeoManager->GetCurrentNodeId()+1; | |
89125dc7 | 1283 | gMCGeom->SetNextRegion(newReg, newLttc); |
efde9b4d | 1284 | flagErr = newReg; |
89125dc7 | 1285 | if (gMCGeom->IsDebugging()) { |
2bc4c610 | 1286 | printf(" out: newReg=%i newLttc=%i\n", newReg, newLttc); |
1287 | printf("<= LKDBWR\n"); | |
1288 | } | |
efde9b4d | 1289 | } |
1290 | ||
1291 | //_____________________________________________________________________________ | |
1292 | void lkfxwr(Double_t &pSx, Double_t &pSy, Double_t &pSz, | |
1293 | Double_t * /*pV*/, const Int_t &oldReg, const Int_t &oldLttc, | |
1294 | Int_t &newReg, Int_t &flagErr, Int_t &newLttc) | |
1295 | { | |
89125dc7 | 1296 | if (gMCGeom->IsDebugging()) { |
2bc4c610 | 1297 | printf("========== Inside LKFXWR (%f, %f, %f)\n",pSx, pSy, pSz); |
1298 | // printf(" in: pV=(%f, %f, %f)\n", pV[0], pV[1], pV[2]); | |
1299 | printf(" in: oldReg=%i oldLttc=%i\n", oldReg, oldLttc); | |
1300 | } | |
efde9b4d | 1301 | TGeoNode *node = gGeoManager->FindNode(pSx, pSy, pSz); |
1302 | if (gGeoManager->IsOutside()) { | |
89125dc7 | 1303 | newReg = gMCGeom->NofVolumes()+1; |
05265ca9 | 1304 | // newLttc = gGeoManager->GetCurrentNodeId(); |
1305 | newLttc = 999999999; | |
89125dc7 | 1306 | if (gMCGeom->IsDebugging()) { |
2bc4c610 | 1307 | printf("OUTSIDE\n"); |
1308 | printf(" out: newReg=%i newLttc=%i\n", newReg, newLttc); | |
1309 | printf("<= LKMGWR\n"); | |
1310 | } | |
05265ca9 | 1311 | flagErr = newReg; |
1312 | return; | |
efde9b4d | 1313 | } |
1314 | newReg = node->GetVolume()->GetNumber(); | |
1315 | newLttc = gGeoManager->GetCurrentNodeId()+1; | |
89125dc7 | 1316 | gMCGeom->SetNextRegion(newReg, newLttc); |
efde9b4d | 1317 | flagErr = newReg; |
89125dc7 | 1318 | if (gMCGeom->IsDebugging()) { |
2bc4c610 | 1319 | printf(" out: newReg=%i newLttc=%i\n", newReg, newLttc); |
1320 | printf("<= LKFXWR\n"); | |
1321 | } | |
efde9b4d | 1322 | } |
1323 | ||
1324 | //_____________________________________________________________________________ | |
1325 | void lkmgwr(Double_t &pSx, Double_t &pSy, Double_t &pSz, | |
1326 | Double_t * /*pV*/, const Int_t &oldReg, const Int_t &oldLttc, | |
1327 | Int_t &flagErr, Int_t &newReg, Int_t &newLttc) | |
1328 | { | |
89125dc7 | 1329 | if (gMCGeom->IsDebugging()) { |
2bc4c610 | 1330 | printf("========== Inside LKMGWR (%f, %f, %f)\n",pSx, pSy, pSz); |
1331 | // printf(" in: pV=(%f, %f, %f)\n", pV[0], pV[1], pV[2]); | |
1332 | printf(" in: oldReg=%i oldLttc=%i\n", oldReg, oldLttc); | |
1333 | } | |
efde9b4d | 1334 | TGeoNode *node = gGeoManager->FindNode(pSx, pSy, pSz); |
1335 | if (gGeoManager->IsOutside()) { | |
89125dc7 | 1336 | newReg = gMCGeom->NofVolumes()+1; |
05265ca9 | 1337 | // newLttc = gGeoManager->GetCurrentNodeId(); |
1338 | newLttc = 999999999; | |
89125dc7 | 1339 | if (gMCGeom->IsDebugging()) { |
2bc4c610 | 1340 | printf("OUTSIDE\n"); |
1341 | printf(" out: newReg=%i newLttc=%i\n", newReg, newLttc); | |
1342 | printf("<= LKMGWR\n"); | |
1343 | } | |
05265ca9 | 1344 | flagErr = newReg; |
1345 | return; | |
efde9b4d | 1346 | } |
1347 | newReg = node->GetVolume()->GetNumber(); | |
1348 | newLttc = gGeoManager->GetCurrentNodeId()+1; | |
89125dc7 | 1349 | gMCGeom->SetNextRegion(newReg, newLttc); |
efde9b4d | 1350 | flagErr = newReg; |
89125dc7 | 1351 | if (gMCGeom->IsDebugging()) { |
2bc4c610 | 1352 | printf(" out: newReg=%i newLttc=%i\n", newReg, newLttc); |
1353 | printf("<= LKMGWR\n"); | |
1354 | } | |
efde9b4d | 1355 | } |
1356 | ||
1357 | //_____________________________________________________________________________ | |
1358 | void lkwr(Double_t &pSx, Double_t &pSy, Double_t &pSz, | |
1359 | Double_t * /*pV*/, const Int_t &oldReg, const Int_t &oldLttc, | |
1360 | Int_t &newReg, Int_t &flagErr, Int_t &newLttc) | |
1361 | { | |
89125dc7 | 1362 | if (gMCGeom->IsDebugging()) { |
2bc4c610 | 1363 | printf("========== Inside LKWR (%f, %f, %f)\n",pSx, pSy, pSz); |
1364 | // printf(" in: pV=(%f, %f, %f)\n", pV[0], pV[1], pV[2]); | |
1365 | printf(" in: oldReg=%i oldLttc=%i\n", oldReg, oldLttc); | |
1366 | } | |
efde9b4d | 1367 | TGeoNode *node = gGeoManager->FindNode(pSx, pSy, pSz); |
1368 | if (gGeoManager->IsOutside()) { | |
89125dc7 | 1369 | newReg = gMCGeom->NofVolumes()+1; |
05265ca9 | 1370 | // newLttc = gGeoManager->GetCurrentNodeId(); |
1371 | newLttc = 999999999; | |
89125dc7 | 1372 | if (gMCGeom->IsDebugging()) { |
2bc4c610 | 1373 | printf("OUTSIDE\n"); |
1374 | printf(" out: newReg=%i newLttc=%i\n", newReg, newLttc); | |
1375 | printf("<= LKMGWR\n"); | |
1376 | } | |
05265ca9 | 1377 | flagErr = newReg; |
1378 | return; | |
efde9b4d | 1379 | } |
1380 | newReg = node->GetVolume()->GetNumber(); | |
d11bf3ca | 1381 | newLttc = gGeoManager->GetCurrentNodeId()+1; |
89125dc7 | 1382 | gMCGeom->SetNextRegion(newReg, newLttc); |
efde9b4d | 1383 | flagErr = newReg; |
89125dc7 | 1384 | if (gMCGeom->IsDebugging()) { |
2bc4c610 | 1385 | printf(" out: newReg=%i newLttc=%i in %s\n", newReg, newLttc, gGeoManager->GetPath()); |
1386 | printf("<= LKWR\n"); | |
1387 | } | |
efde9b4d | 1388 | } |
1389 | ||
1390 | //_____________________________________________________________________________ | |
2573ac89 | 1391 | void nrmlwr(Double_t &pSx, Double_t &pSy, Double_t &pSz, |
1392 | Double_t &pVx, Double_t &pVy, Double_t &pVz, | |
1393 | Double_t *norml, const Int_t &oldReg, | |
1394 | const Int_t &newReg, Int_t &flagErr) | |
efde9b4d | 1395 | { |
89125dc7 | 1396 | if (gMCGeom->IsDebugging()) { |
2bc4c610 | 1397 | printf("========== Inside NRMLWR (%g, %g, %g, %g, %g, %g)\n", pSx,pSy,pSz,pVx,pVy,pVz); |
1398 | printf(" oldReg=%i, newReg=%i\n", oldReg,newReg); | |
1399 | } | |
89125dc7 | 1400 | // Int_t curreg = (gGeoManager->IsOutside())?(gMCGeom->NofVolumes()+1):gGeoManager->GetCurrentVolume()->GetNumber(); |
d11bf3ca | 1401 | // Int_t curLttc = gGeoManager->GetCurrentNodeId()+1; |
89125dc7 | 1402 | // if (gMCGeom->IsDebugging()) printf(" curReg=%i, curLttc=%i in: %s\n", curreg, curLttc, gGeoManager->GetPath()); |
d11bf3ca | 1403 | // Bool_t regsame = (curreg==oldReg)?kTRUE:kFALSE; |
2573ac89 | 1404 | gGeoManager->SetCurrentPoint(pSx, pSy, pSz); |
1405 | gGeoManager->SetCurrentDirection(pVx,pVy,pVz); | |
d11bf3ca | 1406 | /* |
2573ac89 | 1407 | if (!regsame) { |
89125dc7 | 1408 | if (gMCGeom->IsDebugging()) printf(" REGIONS DOEN NOT MATCH\n"); |
2573ac89 | 1409 | gGeoManager->FindNode(); |
89125dc7 | 1410 | curreg = (gGeoManager->IsOutside())?(gMCGeom->NofVolumes()+1):gGeoManager->GetCurrentVolume()->GetNumber(); |
2573ac89 | 1411 | curLttc = gGeoManager->GetCurrentNodeId()+1; |
89125dc7 | 1412 | if (gMCGeom->IsDebugging()) printf(" re-initialized point: curReg=%i curLttc=%i curPath=%s\n", curreg, curLttc, gGeoManager->GetPath()); |
2573ac89 | 1413 | } |
d11bf3ca | 1414 | */ |
2573ac89 | 1415 | Double_t *dnorm = gGeoManager->FindNormalFast(); |
1416 | flagErr = 0; | |
1417 | if (!dnorm) { | |
1418 | printf(" ERROR: Cannot compute fast normal\n"); | |
1419 | flagErr = 1; | |
1420 | norml[0] = -pVx; | |
1421 | norml[1] = -pVy; | |
1422 | norml[2] = -pVz; | |
1423 | } | |
1424 | norml[0] = -dnorm[0]; | |
1425 | norml[1] = -dnorm[1]; | |
1426 | norml[2] = -dnorm[2]; | |
89125dc7 | 1427 | if (gMCGeom->IsDebugging()) printf(" normal to boundary: (%g, %g, %g)\n", norml[0], norml[1], norml[2]); |
1428 | // curreg = (gGeoManager->IsOutside())?(gMCGeom->NofVolumes()+1):gGeoManager->GetCurrentVolume()->GetNumber(); | |
d11bf3ca | 1429 | // curLttc = gGeoManager->GetCurrentNodeId()+1; |
89125dc7 | 1430 | if (gMCGeom->IsDebugging()) { |
d11bf3ca | 1431 | // printf(" final location: curReg=%i, curLttc=%i in %s\n", curreg,curLttc,gGeoManager->GetPath()); |
2bc4c610 | 1432 | printf("<= NRMLWR\n"); |
1433 | } | |
efde9b4d | 1434 | } |
1435 | ||
1436 | //_____________________________________________________________________________ | |
1437 | void rgrpwr(const Int_t & /*flukaReg*/, const Int_t & /*ptrLttc*/, Int_t & /*g4Reg*/, | |
1438 | Int_t * /*indMother*/, Int_t * /*repMother*/, Int_t & /*depthFluka*/) | |
1439 | { | |
89125dc7 | 1440 | if (gMCGeom->IsDebugging()) printf("=> Dummy RGRPWR\n"); |
efde9b4d | 1441 | } |
1442 | ||
1443 | //_____________________________________________________________________________ | |
1444 | Int_t isvhwr(const Int_t &check, const Int_t & intHist) | |
1445 | { | |
1446 | // from FLUGG: | |
1447 | // Wrapper for saving current navigation history (fCheck=default) | |
1448 | // and returning its pointer. If fCheck=-1 copy of history pointed | |
1449 | // by intHist is made in NavHistWithCount object, and its pointer | |
1450 | // is returned. fCheck=1 and fCheck=2 cases are only in debugging | |
1451 | // version: an array is created by means of FGeometryInit functions | |
1452 | // (but could be a static int * ptrArray = new int[10000] with | |
1453 | // file scope as well) that stores a flag for deleted/undeleted | |
1454 | // histories and at the end of event is checked to verify that | |
1455 | // all saved history objects have been deleted. | |
1456 | ||
1457 | // For TGeo, just return the current node ID. No copy need to be made. | |
1458 | ||
89125dc7 | 1459 | if (gMCGeom->IsDebugging()) printf("=> Inside ISVHWR\n"); |
efde9b4d | 1460 | if (check<0) return intHist; |
2573ac89 | 1461 | Int_t histInt = gGeoManager->GetCurrentNodeId()+1; |
89125dc7 | 1462 | if (gMCGeom->IsDebugging()) printf("<= ISVHWR: history is: %i in: %s\n", histInt, gGeoManager->GetPath()); |
efde9b4d | 1463 | return histInt; |
1464 | } | |
1465 | ||
1466 | ||
8495a208 | 1467 | |
1468 |