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