]>
Commit | Line | Data |
---|---|---|
8495a208 | 1 | // @(#):$Name$:$Id$ |
2 | // Author: Andrei Gheata 10/07/2003 | |
3 | ||
4 | #include "TObjString.h" | |
d11bf3ca | 5 | #include "TFlukaGeo.h" |
8495a208 | 6 | //#include "TVirtualMCApplication.h" |
7 | #include "TFlukaMCGeometry.h" | |
8 | #include "TGeoManager.h" | |
9 | #include "TGeoVolume.h" | |
10 | ||
11 | #include "TCallf77.h" | |
12 | ||
13 | #ifndef WIN32 | |
14 | # define idnrwr idnrwr_ | |
15 | # define g1wr g1wr_ | |
16 | # define g1rtwr g1rtwr_ | |
17 | # define conhwr conhwr_ | |
18 | # define inihwr inihwr_ | |
19 | # define jomiwr jomiwr_ | |
20 | # define lkdbwr lkdbwr_ | |
21 | # define lkfxwr lkfxwr_ | |
22 | # define lkmgwr lkmgwr_ | |
23 | # define lkwr lkwr_ | |
24 | # define magfld magfld_ | |
25 | # define nrmlwr nrmlwr_ | |
26 | # define rgrpwr rgrpwr_ | |
27 | # define isvhwr isvhwr_ | |
28 | ||
29 | #else | |
30 | ||
31 | # define idnrwr IDNRWR | |
32 | # define g1wr G1WR | |
33 | # define g1rtwr G1RTWR | |
34 | # define conhwr CONHWR | |
35 | # define inihwr INIHWR | |
36 | # define jomiwr JOMIWR | |
37 | # define lkdbwr LKDBWR | |
38 | # define lkfxwr LKFXWR | |
39 | # define lkmgwr LKMGWR | |
40 | # define lkwr LKWR | |
41 | # define magfld MAGFLD | |
42 | # define nrmlwr NRMLWR | |
43 | # define rgrpwr RGRPWR | |
44 | # define isvhwr ISVHWR | |
45 | ||
46 | #endif | |
47 | ||
48 | //____________________________________________________________________________ | |
49 | extern "C" | |
50 | { | |
51 | // | |
52 | // Prototypes for FLUKA navigation methods | |
53 | // | |
54 | Int_t type_of_call idnrwr(const Int_t & /*nreg*/, const Int_t & /*mlat*/); | |
55 | void type_of_call g1wr(Double_t & /*pSx*/, Double_t & /*pSy*/, Double_t & /*pSz*/, | |
56 | Double_t * /*pV*/, Int_t & /*oldReg*/ , const Int_t & /*oldLttc*/, Double_t & /*propStep*/, | |
57 | Int_t & /*nascFlag*/, Double_t & /*retStep*/, Int_t & /*newReg*/, | |
58 | Double_t & /*saf*/, Int_t & /*newLttc*/, Int_t & /*LttcFlag*/, | |
59 | Double_t *s /*Lt*/, Int_t * /*jrLt*/); | |
60 | ||
61 | void type_of_call g1rtwr(); | |
62 | void type_of_call conhwr(Int_t & /*intHist*/, Int_t * /*incrCount*/); | |
63 | void type_of_call inihwr(Int_t & /*intHist*/); | |
64 | void type_of_call jomiwr(const Int_t & /*nge*/, const Int_t & /*lin*/, const Int_t & /*lou*/, | |
65 | Int_t & /*flukaReg*/); | |
66 | void type_of_call lkdbwr(Double_t & /*pSx*/, Double_t & /*pSy*/, Double_t & /*pSz*/, | |
67 | Double_t * /*pV*/, const Int_t & /*oldReg*/, const Int_t & /*oldLttc*/, | |
68 | Int_t & /*newReg*/, Int_t & /*flagErr*/, Int_t & /*newLttc*/); | |
69 | void type_of_call lkfxwr(Double_t & /*pSx*/, Double_t & /*pSy*/, Double_t & /*pSz*/, | |
70 | Double_t * /*pV*/, const Int_t & /*oldReg*/, const Int_t & /*oldLttc*/, | |
71 | Int_t & /*newReg*/, Int_t & /*flagErr*/, Int_t & /*newLttc*/); | |
72 | void type_of_call lkmgwr(Double_t & /*pSx*/, Double_t & /*pSy*/, Double_t & /*pSz*/, | |
73 | Double_t * /*pV*/, const Int_t & /*oldReg*/, const Int_t & /*oldLttc*/, | |
74 | Int_t & /*flagErr*/, Int_t & /*newReg*/, Int_t & /*newLttc*/); | |
75 | void type_of_call lkwr(Double_t & /*pSx*/, Double_t & /*pSy*/, Double_t & /*pSz*/, | |
76 | Double_t * /*pV*/, const Int_t & /*oldReg*/, const Int_t & /*oldLttc*/, | |
77 | Int_t & /*newReg*/, Int_t & /*flagErr*/, Int_t & /*newLttc*/); | |
efde9b4d | 78 | // void type_of_call magfld(const Double_t & /*pX*/, const Double_t & /*pY*/, const Double_t & /*pZ*/, |
79 | // Double_t & /*cosBx*/, Double_t & /*cosBy*/, Double_t & /*cosBz*/, | |
80 | // Double_t & /*Bmag*/, Int_t & /*reg*/, Int_t & /*idiscflag*/); | |
8495a208 | 81 | void type_of_call nrmlwr(Double_t & /*pSx*/, Double_t & /*pSy*/, Double_t & /*pSz*/, |
82 | Double_t & /*pVx*/, Double_t & /*pVy*/, Double_t & /*pVz*/, | |
83 | Double_t * /*norml*/, const Int_t & /*oldReg*/, | |
84 | const Int_t & /*newReg*/, Int_t & /*flagErr*/); | |
85 | void type_of_call rgrpwr(const Int_t & /*flukaReg*/, const Int_t & /*ptrLttc*/, Int_t & /*g4Reg*/, | |
86 | Int_t * /*indMother*/, Int_t * /*repMother*/, Int_t & /*depthFluka*/); | |
87 | Int_t type_of_call isvhwr(const Int_t & /*fCheck*/, const Int_t & /*intHist*/); | |
88 | }; | |
89 | ||
90 | // TFluka global pointer | |
91 | TFluka *fluka = 0; | |
92 | TFlukaMCGeometry *mcgeom = 0; | |
d11bf3ca | 93 | Int_t kNstep = 0; |
8495a208 | 94 | |
95 | ClassImp(TFlukaMCGeometry) | |
96 | ||
97 | TFlukaMCGeometry* TFlukaMCGeometry::fgInstance=0; | |
98 | ||
99 | //_____________________________________________________________________________ | |
100 | TFlukaMCGeometry::TFlukaMCGeometry(const char *name, const char *title) | |
101 | : TVirtualMCGeometry(name, title) | |
102 | { | |
103 | // | |
104 | // Standard constructor | |
105 | // | |
2bc4c610 | 106 | fDebug = kFALSE; |
efde9b4d | 107 | fLastMaterial = 0; |
b1536e91 | 108 | fNextRegion = 0; |
109 | fNextLattice = 0; | |
110 | fRegionList = 0; | |
8495a208 | 111 | fluka = (TFluka*)gMC; |
112 | mcgeom = this; | |
d11bf3ca | 113 | kNstep = 0; |
8495a208 | 114 | } |
115 | ||
116 | //_____________________________________________________________________________ | |
117 | TFlukaMCGeometry::TFlukaMCGeometry() | |
118 | : TVirtualMCGeometry() | |
119 | { | |
120 | // | |
121 | // Default constructor | |
122 | // | |
2bc4c610 | 123 | fDebug = kFALSE; |
efde9b4d | 124 | fLastMaterial = 0; |
b1536e91 | 125 | fNextRegion = 0; |
126 | fNextLattice = 0; | |
127 | fRegionList = 0; | |
8495a208 | 128 | fluka = (TFluka*)gMC; |
129 | mcgeom = this; | |
d11bf3ca | 130 | kNstep = 0; |
8495a208 | 131 | } |
132 | ||
133 | //_____________________________________________________________________________ | |
134 | TFlukaMCGeometry::~TFlukaMCGeometry() | |
135 | { | |
136 | // | |
137 | // Destructor | |
138 | // | |
139 | fgInstance=0; | |
b1536e91 | 140 | if (fRegionList) delete [] fRegionList; |
8495a208 | 141 | if (gGeoManager) delete gGeoManager; |
142 | } | |
143 | ||
144 | // | |
145 | // private methods | |
146 | // | |
147 | //_____________________________________________________________________________ | |
148 | TFlukaMCGeometry::TFlukaMCGeometry(const TFlukaMCGeometry &) | |
149 | : TVirtualMCGeometry() | |
150 | { | |
151 | // | |
152 | // Copy constructor | |
153 | // | |
154 | } | |
155 | ||
156 | //_____________________________________________________________________________ | |
157 | Double_t* TFlukaMCGeometry::CreateDoubleArray(Float_t* array, Int_t size) const | |
158 | { | |
159 | // Converts Float_t* array to Double_t*, | |
160 | // !! The new array has to be deleted by user. | |
161 | // --- | |
162 | ||
163 | Double_t* doubleArray; | |
164 | if (size>0) { | |
165 | doubleArray = new Double_t[size]; | |
166 | for (Int_t i=0; i<size; i++) doubleArray[i] = array[i]; | |
167 | } | |
168 | else { | |
169 | //doubleArray = 0; | |
170 | doubleArray = new Double_t[1]; | |
171 | } | |
172 | return doubleArray; | |
173 | } | |
174 | // | |
175 | // public methods | |
176 | //_____________________________________________________________________________ | |
177 | void TFlukaMCGeometry::Gfmate(Int_t imat, char *name, Float_t &a, Float_t &z, | |
178 | Float_t &dens, Float_t &radl, Float_t &absl, | |
179 | Float_t* /*ubuf*/, Int_t& /*nbuf*/) | |
180 | { | |
2bc4c610 | 181 | if (fDebug) printf("Gfmate %i\n", imat); |
8495a208 | 182 | TGeoMaterial *mat; |
183 | TIter next (gGeoManager->GetListOfMaterials()); | |
184 | while ((mat = (TGeoMaterial*)next())) { | |
185 | if (mat->GetUniqueID() == (UInt_t)imat) break; | |
186 | } | |
187 | if (!mat) { | |
188 | Error("Gfmate", "no material with index %i found", imat); | |
189 | return; | |
190 | } | |
191 | sprintf(name, "%s", mat->GetName()); | |
192 | a = mat->GetA(); | |
193 | z = mat->GetZ(); | |
194 | dens = mat->GetDensity(); | |
195 | radl = mat->GetRadLen(); | |
196 | absl = mat->GetIntLen(); | |
2bc4c610 | 197 | if (fDebug) printf(" ->material found : %s a=%g, z=%g, dens=%g, radl=%g, absl=%g\n", name, a,z,dens,radl,absl); |
8495a208 | 198 | } |
199 | ||
200 | //_____________________________________________________________________________ | |
201 | void TFlukaMCGeometry::Gfmate(Int_t imat, char *name, Double_t &a, Double_t &z, | |
202 | Double_t &dens, Double_t &radl, Double_t &absl, | |
203 | Double_t* /*ubuf*/, Int_t& /*nbuf*/) | |
204 | { | |
2bc4c610 | 205 | if (fDebug) printf("Gfmate %i\n", imat); |
206 | TGeoMaterial *mat; | |
8495a208 | 207 | TIter next (gGeoManager->GetListOfMaterials()); |
208 | while ((mat = (TGeoMaterial*)next())) { | |
209 | if (mat->GetUniqueID() == (UInt_t)imat) break; | |
210 | } | |
211 | if (!mat) { | |
212 | Error("Gfmate", "no material with index %i found", imat); | |
213 | return; | |
214 | } | |
215 | sprintf(name, "%s", mat->GetName()); | |
216 | a = mat->GetA(); | |
217 | z = mat->GetZ(); | |
218 | dens = mat->GetDensity(); | |
219 | radl = mat->GetRadLen(); | |
220 | absl = mat->GetIntLen(); | |
2bc4c610 | 221 | if (fDebug) printf(" ->material found : %s a=%g, z=%g, dens=%g, radl=%g, absl=%g\n", name, a,z,dens,radl,absl); |
8495a208 | 222 | } |
223 | ||
224 | //_____________________________________________________________________________ | |
225 | void TFlukaMCGeometry::Material(Int_t& kmat, const char* name, Double_t a, Double_t z, | |
226 | Double_t dens, Double_t radl, Double_t absl, Float_t* buf, | |
227 | Int_t nwbuf) | |
228 | { | |
229 | // | |
230 | // Defines a Material | |
231 | // | |
232 | // kmat number assigned to the material | |
233 | // name material name | |
234 | // a atomic mass in au | |
235 | // z atomic number | |
236 | // dens density in g/cm3 | |
237 | // absl absorbtion length in cm | |
238 | // if >=0 it is ignored and the program | |
239 | // calculates it, if <0. -absl is taken | |
240 | // radl radiation length in cm | |
241 | // if >=0 it is ignored and the program | |
242 | // calculates it, if <0. -radl is taken | |
243 | // buf pointer to an array of user words | |
244 | // nbuf number of user words | |
245 | // | |
246 | ||
247 | Double_t* dbuf = CreateDoubleArray(buf, nwbuf); | |
248 | Material(kmat, name, a, z, dens, radl, absl, dbuf, nwbuf); | |
249 | delete [] dbuf; | |
250 | } | |
251 | ||
252 | //_____________________________________________________________________________ | |
253 | void TFlukaMCGeometry::Material(Int_t& kmat, const char* name, Double_t a, Double_t z, | |
254 | Double_t dens, Double_t radl, Double_t absl, Double_t* /*buf*/, | |
255 | Int_t /*nwbuf*/) | |
256 | { | |
257 | // | |
258 | // Defines a Material | |
259 | // | |
260 | // kmat number assigned to the material | |
261 | // name material name | |
262 | // a atomic mass in au | |
263 | // z atomic number | |
264 | // dens density in g/cm3 | |
265 | // absl absorbtion length in cm | |
266 | // if >=0 it is ignored and the program | |
267 | // calculates it, if <0. -absl is taken | |
268 | // radl radiation length in cm | |
269 | // if >=0 it is ignored and the program | |
270 | // calculates it, if <0. -radl is taken | |
271 | // buf pointer to an array of user words | |
272 | // nbuf number of user words | |
273 | // | |
274 | ||
275 | kmat = gGeoManager->GetListOfMaterials()->GetSize(); | |
276 | gGeoManager->Material(name, a, z, dens, kmat, radl, absl); | |
2bc4c610 | 277 | if (fDebug) printf("Material %s: kmat=%i, a=%g, z=%g, dens=%g\n", name, kmat, a, z, dens); |
8495a208 | 278 | } |
279 | ||
280 | //_____________________________________________________________________________ | |
281 | void TFlukaMCGeometry::Mixture(Int_t& kmat, const char* name, Float_t* a, Float_t* z, | |
282 | Double_t dens, Int_t nlmat, Float_t* wmat) | |
283 | { | |
284 | // | |
285 | // Defines mixture OR COMPOUND IMAT as composed by | |
286 | // THE BASIC NLMAT materials defined by arrays A,Z and WMAT | |
287 | // | |
288 | // If NLMAT > 0 then wmat contains the proportion by | |
289 | // weights of each basic material in the mixture. | |
290 | // | |
291 | // If nlmat < 0 then WMAT contains the number of atoms | |
292 | // of a given kind into the molecule of the COMPOUND | |
293 | // In this case, WMAT in output is changed to relative | |
294 | // weigths. | |
295 | // | |
296 | ||
297 | Double_t* da = CreateDoubleArray(a, TMath::Abs(nlmat)); | |
298 | Double_t* dz = CreateDoubleArray(z, TMath::Abs(nlmat)); | |
299 | Double_t* dwmat = CreateDoubleArray(wmat, TMath::Abs(nlmat)); | |
300 | ||
301 | Mixture(kmat, name, da, dz, dens, nlmat, dwmat); | |
302 | for (Int_t i=0; i<nlmat; i++) { | |
303 | a[i] = da[i]; z[i] = dz[i]; wmat[i] = dwmat[i]; | |
304 | } | |
305 | ||
306 | delete [] da; | |
307 | delete [] dz; | |
308 | delete [] dwmat; | |
309 | } | |
310 | ||
311 | //_____________________________________________________________________________ | |
312 | void TFlukaMCGeometry::Mixture(Int_t& kmat, const char* name, Double_t* a, Double_t* z, | |
313 | Double_t dens, Int_t nlmat, Double_t* wmat) | |
314 | { | |
315 | // | |
316 | // Defines mixture OR COMPOUND IMAT as composed by | |
317 | // THE BASIC NLMAT materials defined by arrays A,Z and WMAT | |
318 | // | |
319 | // If NLMAT > 0 then wmat contains the proportion by | |
320 | // weights of each basic material in the mixture. | |
321 | // | |
322 | // If nlmat < 0 then WMAT contains the number of atoms | |
323 | // of a given kind into the molecule of the COMPOUND | |
324 | // In this case, WMAT in output is changed to relative | |
325 | // weigths. | |
326 | // | |
327 | ||
328 | if (nlmat < 0) { | |
329 | nlmat = - nlmat; | |
330 | Double_t amol = 0; | |
331 | Int_t i; | |
332 | for (i=0;i<nlmat;i++) { | |
333 | amol += a[i]*wmat[i]; | |
334 | } | |
335 | for (i=0;i<nlmat;i++) { | |
336 | wmat[i] *= a[i]/amol; | |
337 | } | |
338 | } | |
339 | kmat = gGeoManager->GetListOfMaterials()->GetSize(); | |
2bc4c610 | 340 | if (fDebug) { |
341 | printf("Mixture %s with %i elem: kmat=%i, dens=%g\n", name, nlmat, kmat, dens); | |
342 | for (Int_t j=0; j<nlmat; j++) printf(" Elem %i: z=%g a=%g w=%g\n",j,z[j],a[j],wmat[j]); | |
343 | } | |
8495a208 | 344 | gGeoManager->Mixture(name, a, z, dens, nlmat, wmat, kmat); |
345 | } | |
346 | //_____________________________________________________________________________ | |
347 | Int_t TFlukaMCGeometry::GetMedium() const | |
348 | { | |
349 | // Get current medium number | |
350 | Int_t imed = 0; | |
351 | TGeoNode *node = gGeoManager->GetCurrentNode(); | |
352 | if (!node) imed = gGeoManager->GetTopNode()->GetVolume()->GetMedium()->GetId(); | |
353 | else imed = node->GetVolume()->GetMedium()->GetId(); | |
2bc4c610 | 354 | if (fDebug) printf("GetMedium=%i\n", imed); |
8495a208 | 355 | return imed; |
356 | } | |
357 | ||
0bc73b6c | 358 | //_____________________________________________________________________________ |
359 | Int_t TFlukaMCGeometry::GetFlukaMaterial(Int_t imed) const | |
360 | { | |
361 | // Returns FLUKA material index for medium IMED | |
362 | TGeoMedium *med = (TGeoMedium*)gGeoManager->GetListOfMedia()->At(imed-1); | |
363 | if (!med) { | |
364 | Error("GetFlukaMaterial", "MEDIUM %i nor found", imed); | |
365 | return -1; | |
366 | } | |
367 | Int_t imatfl = med->GetMaterial()->GetIndex(); | |
368 | return imatfl; | |
369 | } | |
370 | ||
b1536e91 | 371 | //_____________________________________________________________________________ |
372 | Int_t *TFlukaMCGeometry::GetRegionList(Int_t imed, Int_t &nreg) | |
373 | { | |
374 | // Get an ordered list of regions matching a given medium number | |
375 | nreg = 0; | |
376 | if (!fRegionList) fRegionList = new Int_t[NofVolumes()+1]; | |
377 | TIter next(gGeoManager->GetListOfUVolumes()); | |
378 | TGeoVolume *vol; | |
379 | Int_t imedium, ireg; | |
380 | while ((vol = (TGeoVolume*)next())) { | |
381 | imedium = vol->GetMedium()->GetId(); | |
382 | if (imedium == imed) { | |
383 | ireg = vol->GetNumber(); | |
384 | fRegionList[nreg++] = ireg; | |
385 | } | |
386 | } | |
387 | return fRegionList; | |
388 | } | |
389 | ||
390 | //_____________________________________________________________________________ | |
391 | Int_t *TFlukaMCGeometry::GetMaterialList(Int_t imat, Int_t &nreg) | |
392 | { | |
393 | // Get an ordered list of regions matching a given medium number | |
394 | nreg = 0; | |
395 | if (!fRegionList) fRegionList = new Int_t[NofVolumes()+1]; | |
396 | TIter next(gGeoManager->GetListOfUVolumes()); | |
397 | TGeoVolume *vol; | |
398 | Int_t imaterial, ireg; | |
399 | while ((vol = (TGeoVolume*)next())) { | |
400 | imaterial = vol->GetMedium()->GetMaterial()->GetIndex(); | |
401 | if (imaterial == imat) { | |
402 | ireg = vol->GetNumber(); | |
403 | fRegionList[nreg++] = ireg; | |
404 | } | |
405 | } | |
406 | return fRegionList; | |
407 | } | |
8495a208 | 408 | //_____________________________________________________________________________ |
409 | void TFlukaMCGeometry::Medium(Int_t& kmed, const char* name, Int_t nmat, Int_t isvol, | |
410 | Int_t ifield, Double_t fieldm, Double_t tmaxfd, | |
411 | Double_t stemax, Double_t deemax, Double_t epsil, | |
412 | Double_t stmin, Float_t* ubuf, Int_t nbuf) | |
413 | { | |
414 | // | |
415 | // kmed tracking medium number assigned | |
416 | // name tracking medium name | |
417 | // nmat material number | |
418 | // isvol sensitive volume flag | |
419 | // ifield magnetic field | |
420 | // fieldm max. field value (kilogauss) | |
421 | // tmaxfd max. angle due to field (deg/step) | |
422 | // stemax max. step allowed | |
423 | // deemax max. fraction of energy lost in a step | |
424 | // epsil tracking precision (cm) | |
425 | // stmin min. step due to continuous processes (cm) | |
426 | // | |
427 | // ifield = 0 if no magnetic field; ifield = -1 if user decision in guswim; | |
428 | // ifield = 1 if tracking performed with g3rkuta; ifield = 2 if tracking | |
429 | // performed with g3helix; ifield = 3 if tracking performed with g3helx3. | |
430 | // | |
431 | ||
432 | //printf("Creating mediuma: %s, numed=%d, nmat=%d\n",name,kmed,nmat); | |
433 | Double_t* dubuf = CreateDoubleArray(ubuf, nbuf); | |
434 | Medium(kmed, name, nmat, isvol, ifield, fieldm, tmaxfd, stemax, deemax, epsil, | |
435 | stmin, dubuf, nbuf); | |
436 | delete [] dubuf; | |
437 | } | |
438 | ||
439 | //_____________________________________________________________________________ | |
440 | void TFlukaMCGeometry::Medium(Int_t& kmed, const char* name, Int_t nmat, Int_t isvol, | |
441 | Int_t ifield, Double_t fieldm, Double_t tmaxfd, | |
442 | Double_t stemax, Double_t deemax, Double_t epsil, | |
443 | Double_t stmin, Double_t* /*ubuf*/, Int_t /*nbuf*/) | |
444 | { | |
445 | // | |
446 | // kmed tracking medium number assigned | |
447 | // name tracking medium name | |
448 | // nmat material number | |
449 | // isvol sensitive volume flag | |
450 | // ifield magnetic field | |
451 | // fieldm max. field value (kilogauss) | |
452 | // tmaxfd max. angle due to field (deg/step) | |
453 | // stemax max. step allowed | |
454 | // deemax max. fraction of energy lost in a step | |
455 | // epsil tracking precision (cm) | |
456 | // stmin min. step due to continuos processes (cm) | |
457 | // | |
458 | // ifield = 0 if no magnetic field; ifield = -1 if user decision in guswim; | |
459 | // ifield = 1 if tracking performed with g3rkuta; ifield = 2 if tracking | |
460 | // performed with g3helix; ifield = 3 if tracking performed with g3helx3. | |
461 | // | |
462 | ||
0bc73b6c | 463 | kmed = gGeoManager->GetListOfMedia()->GetSize()+1; |
8495a208 | 464 | gGeoManager->Medium(name,kmed,nmat, isvol, ifield, fieldm, tmaxfd, stemax,deemax, epsil, stmin); |
2bc4c610 | 465 | if (fDebug) printf("Medium %s: kmed=%i, nmat=%i, isvol=%i\n", name, kmed, nmat,isvol); |
8495a208 | 466 | } |
467 | ||
468 | //_____________________________________________________________________________ | |
469 | void TFlukaMCGeometry::Matrix(Int_t& krot, Double_t thex, Double_t phix, Double_t they, | |
470 | Double_t phiy, Double_t thez, Double_t phiz) | |
471 | { | |
472 | // | |
473 | // krot rotation matrix number assigned | |
474 | // theta1 polar angle for axis i | |
475 | // phi1 azimuthal angle for axis i | |
476 | // theta2 polar angle for axis ii | |
477 | // phi2 azimuthal angle for axis ii | |
478 | // theta3 polar angle for axis iii | |
479 | // phi3 azimuthal angle for axis iii | |
480 | // | |
481 | // it defines the rotation matrix number irot. | |
482 | // | |
483 | ||
484 | krot = gGeoManager->GetListOfMatrices()->GetEntriesFast(); | |
485 | gGeoManager->Matrix(krot, thex, phix, they, phiy, thez, phiz); | |
2bc4c610 | 486 | if (fDebug) printf("Rotation %i defined\n", krot); |
8495a208 | 487 | } |
488 | ||
489 | //_____________________________________________________________________________ | |
490 | Int_t TFlukaMCGeometry::Gsvolu(const char *name, const char *shape, Int_t nmed, | |
491 | Float_t *upar, Int_t npar) | |
492 | { | |
493 | // | |
494 | // NAME Volume name | |
495 | // SHAPE Volume type | |
496 | // NUMED Tracking medium number | |
497 | // NPAR Number of shape parameters | |
498 | // UPAR Vector containing shape parameters | |
499 | // | |
500 | // It creates a new volume in the JVOLUM data structure. | |
501 | // | |
502 | ||
503 | Double_t* dupar = CreateDoubleArray(upar, npar); | |
504 | Int_t id = Gsvolu(name, shape, nmed, dupar, npar); | |
505 | delete [] dupar; | |
506 | return id; | |
507 | } | |
508 | ||
509 | //_____________________________________________________________________________ | |
510 | Int_t TFlukaMCGeometry::Gsvolu(const char *name, const char *shape, Int_t nmed, | |
511 | Double_t *upar, Int_t npar) | |
512 | { | |
513 | // | |
514 | // NAME Volume name | |
515 | // SHAPE Volume type | |
516 | // NUMED Tracking medium number | |
517 | // NPAR Number of shape parameters | |
518 | // UPAR Vector containing shape parameters | |
519 | // | |
520 | // It creates a new volume in the JVOLUM data structure. | |
521 | // | |
522 | char vname[5]; | |
523 | Vname(name,vname); | |
524 | char vshape[5]; | |
525 | Vname(shape,vshape); | |
526 | ||
527 | TGeoVolume* vol = gGeoManager->Volume(vname, shape, nmed, upar, npar); | |
2bc4c610 | 528 | if (fDebug) printf("Volume %s: id=%i shape=%s, nmed=%i\n", vname, vol->GetNumber(), shape, nmed); |
8495a208 | 529 | return vol->GetNumber(); |
530 | } | |
531 | ||
532 | //_____________________________________________________________________________ | |
533 | void TFlukaMCGeometry::Gsdvn(const char *name, const char *mother, Int_t ndiv, | |
534 | Int_t iaxis) | |
535 | { | |
536 | // | |
537 | // Create a new volume by dividing an existing one | |
538 | // | |
539 | // NAME Volume name | |
540 | // MOTHER Mother volume name | |
541 | // NDIV Number of divisions | |
542 | // IAXIS Axis value | |
543 | // | |
544 | // X,Y,Z of CAXIS will be translated to 1,2,3 for IAXIS. | |
545 | // It divides a previously defined volume. | |
546 | // | |
547 | char vname[5]; | |
548 | Vname(name,vname); | |
549 | char vmother[5]; | |
550 | Vname(mother,vmother); | |
551 | ||
552 | gGeoManager->Division(vname, vmother, iaxis, ndiv, 0, 0, 0, "n"); | |
2bc4c610 | 553 | if (fDebug) printf("Division %s: mother=%s iaxis=%i ndiv=%i\n", vname, vmother, iaxis, ndiv); |
8495a208 | 554 | } |
555 | ||
556 | //_____________________________________________________________________________ | |
557 | void TFlukaMCGeometry::Gsdvn2(const char *name, const char *mother, Int_t ndiv, | |
558 | Int_t iaxis, Double_t c0i, Int_t numed) | |
559 | { | |
560 | // | |
561 | // Create a new volume by dividing an existing one | |
562 | // | |
563 | // Divides mother into ndiv divisions called name | |
564 | // along axis iaxis starting at coordinate value c0. | |
565 | // the new volume created will be medium number numed. | |
566 | // | |
567 | char vname[5]; | |
568 | Vname(name,vname); | |
569 | char vmother[5]; | |
570 | Vname(mother,vmother); | |
571 | ||
572 | gGeoManager->Division(vname, vmother, iaxis, ndiv, c0i, 0, numed, "nx"); | |
573 | } | |
574 | //_____________________________________________________________________________ | |
575 | void TFlukaMCGeometry::Gsdvt(const char *name, const char *mother, Double_t step, | |
576 | Int_t iaxis, Int_t numed, Int_t /*ndvmx*/) | |
577 | { | |
578 | // | |
579 | // Create a new volume by dividing an existing one | |
580 | // | |
581 | // Divides MOTHER into divisions called NAME along | |
582 | // axis IAXIS in steps of STEP. If not exactly divisible | |
583 | // will make as many as possible and will centre them | |
584 | // with respect to the mother. Divisions will have medium | |
585 | // number NUMED. If NUMED is 0, NUMED of MOTHER is taken. | |
586 | // NDVMX is the expected maximum number of divisions | |
587 | // (If 0, no protection tests are performed) | |
588 | // | |
589 | char vname[5]; | |
590 | Vname(name,vname); | |
591 | char vmother[5]; | |
592 | Vname(mother,vmother); | |
593 | ||
594 | gGeoManager->Division(vname, vmother, iaxis, 0, 0, step, numed, "s"); | |
595 | } | |
596 | ||
597 | //_____________________________________________________________________________ | |
598 | void TFlukaMCGeometry::Gsdvt2(const char *name, const char *mother, Double_t step, | |
599 | Int_t iaxis, Double_t c0, Int_t numed, Int_t /*ndvmx*/) | |
600 | { | |
601 | // | |
602 | // Create a new volume by dividing an existing one | |
603 | // | |
604 | // Divides MOTHER into divisions called NAME along | |
605 | // axis IAXIS starting at coordinate value C0 with step | |
606 | // size STEP. | |
607 | // The new volume created will have medium number NUMED. | |
608 | // If NUMED is 0, NUMED of mother is taken. | |
609 | // NDVMX is the expected maximum number of divisions | |
610 | // (If 0, no protection tests are performed) | |
611 | // | |
612 | char vname[5]; | |
613 | Vname(name,vname); | |
614 | char vmother[5]; | |
615 | Vname(mother,vmother); | |
616 | ||
617 | gGeoManager->Division(vname, vmother, iaxis, 0, c0, step, numed, "sx"); | |
618 | } | |
619 | ||
620 | //_____________________________________________________________________________ | |
621 | void TFlukaMCGeometry::Gsord(const char * /*name*/, Int_t /*iax*/) | |
622 | { | |
623 | // | |
624 | // Flags volume CHNAME whose contents will have to be ordered | |
625 | // along axis IAX, by setting the search flag to -IAX | |
626 | // IAX = 1 X axis | |
627 | // IAX = 2 Y axis | |
628 | // IAX = 3 Z axis | |
629 | // IAX = 4 Rxy (static ordering only -> GTMEDI) | |
630 | // IAX = 14 Rxy (also dynamic ordering -> GTNEXT) | |
631 | // IAX = 5 Rxyz (static ordering only -> GTMEDI) | |
632 | // IAX = 15 Rxyz (also dynamic ordering -> GTNEXT) | |
633 | // IAX = 6 PHI (PHI=0 => X axis) | |
634 | // IAX = 7 THETA (THETA=0 => Z axis) | |
635 | // | |
636 | ||
637 | // TBC - keep this function | |
638 | // nothing to be done for TGeo //xx | |
639 | } | |
640 | ||
641 | //_____________________________________________________________________________ | |
642 | void TFlukaMCGeometry::Gspos(const char *name, Int_t nr, const char *mother, Double_t x, | |
643 | Double_t y, Double_t z, Int_t irot, const char *konly) | |
644 | { | |
645 | // | |
646 | // Position a volume into an existing one | |
647 | // | |
648 | // NAME Volume name | |
649 | // NUMBER Copy number of the volume | |
650 | // MOTHER Mother volume name | |
651 | // X X coord. of the volume in mother ref. sys. | |
652 | // Y Y coord. of the volume in mother ref. sys. | |
653 | // Z Z coord. of the volume in mother ref. sys. | |
654 | // IROT Rotation matrix number w.r.t. mother ref. sys. | |
655 | // ONLY ONLY/MANY flag | |
656 | // | |
657 | // It positions a previously defined volume in the mother. | |
658 | // | |
659 | ||
660 | TString only = konly; | |
661 | only.ToLower(); | |
662 | Bool_t isOnly = kFALSE; | |
663 | if (only.Contains("only")) isOnly = kTRUE; | |
664 | char vname[5]; | |
665 | Vname(name,vname); | |
666 | char vmother[5]; | |
667 | Vname(mother,vmother); | |
668 | ||
669 | Double_t *upar=0; | |
670 | gGeoManager->Node(vname, nr, vmother, x, y, z, irot, isOnly, upar); | |
2bc4c610 | 671 | if (fDebug) printf("Adding daughter %s to %s: cpy=%i irot=%i only=%s\n", vname,vmother,nr,irot,only.Data()); |
8495a208 | 672 | } |
673 | ||
674 | //_____________________________________________________________________________ | |
675 | void TFlukaMCGeometry::Gsposp(const char *name, Int_t nr, const char *mother, | |
676 | Double_t x, Double_t y, Double_t z, Int_t irot, | |
677 | const char *konly, Float_t *upar, Int_t np ) | |
678 | { | |
679 | // | |
680 | // Place a copy of generic volume NAME with user number | |
681 | // NR inside MOTHER, with its parameters UPAR(1..NP) | |
682 | // | |
683 | ||
684 | Double_t* dupar = CreateDoubleArray(upar, np); | |
685 | Gsposp(name, nr, mother, x, y, z, irot, konly, dupar, np); | |
686 | delete [] dupar; | |
687 | } | |
688 | ||
689 | //_____________________________________________________________________________ | |
690 | void TFlukaMCGeometry::Gsposp(const char *name, Int_t nr, const char *mother, | |
691 | Double_t x, Double_t y, Double_t z, Int_t irot, | |
692 | const char *konly, Double_t *upar, Int_t np ) | |
693 | { | |
694 | // | |
695 | // Place a copy of generic volume NAME with user number | |
696 | // NR inside MOTHER, with its parameters UPAR(1..NP) | |
697 | // | |
698 | ||
699 | TString only = konly; | |
700 | only.ToLower(); | |
701 | Bool_t isOnly = kFALSE; | |
702 | if (only.Contains("only")) isOnly = kTRUE; | |
703 | char vname[5]; | |
704 | Vname(name,vname); | |
705 | char vmother[5]; | |
706 | Vname(mother,vmother); | |
707 | ||
708 | gGeoManager->Node(vname,nr,vmother, x,y,z,irot,isOnly,upar,np); | |
2bc4c610 | 709 | if (fDebug) printf("Adding daughter(s) %s to %s: cpy=%i irot=%i only=%s\n", vname,vmother,nr,irot,only.Data()); |
8495a208 | 710 | } |
711 | ||
712 | //_____________________________________________________________________________ | |
713 | Int_t TFlukaMCGeometry::VolId(const Text_t *name) const | |
714 | { | |
715 | // | |
716 | // Return the unique numeric identifier for volume name | |
717 | // | |
718 | ||
719 | Int_t uid = gGeoManager->GetUID(name); | |
720 | if (uid<0) { | |
721 | printf("VolId: Volume %s not found\n",name); | |
722 | return 0; | |
723 | } | |
2bc4c610 | 724 | if (fDebug) printf("VolId for %s: %i\n", name, uid); |
8495a208 | 725 | return uid; |
726 | } | |
727 | ||
728 | //_____________________________________________________________________________ | |
729 | const char* TFlukaMCGeometry::VolName(Int_t id) const | |
730 | { | |
731 | // | |
732 | // Return the volume name given the volume identifier | |
733 | // | |
734 | TGeoVolume *volume = gGeoManager->GetVolume(id); | |
735 | if (!volume) { | |
736 | Error("VolName","volume with id=%d does not exist",id); | |
737 | return "NULL"; | |
738 | } | |
2bc4c610 | 739 | if (fDebug) printf("VolName for id=%i: %s\n", id, volume->GetName()); |
8495a208 | 740 | return volume->GetName(); |
741 | } | |
742 | ||
743 | //_____________________________________________________________________________ | |
744 | Int_t TFlukaMCGeometry::NofVolumes() const | |
745 | { | |
746 | // | |
747 | // Return total number of volumes in the geometry | |
748 | // | |
749 | ||
750 | return gGeoManager->GetListOfUVolumes()->GetEntriesFast()-1; | |
751 | } | |
752 | ||
753 | //_____________________________________________________________________________ | |
754 | Int_t TFlukaMCGeometry::VolId2Mate(Int_t id) const | |
755 | { | |
756 | // | |
757 | // Return material number for a given volume id | |
758 | // | |
759 | TGeoVolume *volume = gGeoManager->GetVolume(id); | |
760 | if (!volume) { | |
761 | Error("VolId2Mate","volume with id=%d does not exist",id); | |
762 | return 0; | |
763 | } | |
764 | TGeoMedium *med = volume->GetMedium(); | |
765 | if (!med) return 0; | |
2bc4c610 | 766 | if (fDebug) printf("VolId2Mate id=%i: idmed=%i\n", id, med->GetId()); |
8495a208 | 767 | return med->GetId(); |
768 | } | |
769 | ||
770 | //_____________________________________________________________________________ | |
771 | Int_t TFlukaMCGeometry::CurrentVolID(Int_t& copyNo) const | |
772 | { | |
773 | // Returns the current volume ID and copy number | |
774 | if (gGeoManager->IsOutside()) return 0; | |
775 | TGeoNode *node = gGeoManager->GetCurrentNode(); | |
776 | copyNo = node->GetNumber(); | |
777 | Int_t id = node->GetVolume()->GetNumber(); | |
2bc4c610 | 778 | if (fDebug) printf("CurrentVolId(cpy=%i) = %i\n", copyNo, id); |
8495a208 | 779 | return id; |
780 | } | |
781 | ||
782 | //_____________________________________________________________________________ | |
783 | Int_t TFlukaMCGeometry::CurrentVolOffID(Int_t off, Int_t& copyNo) const | |
784 | { | |
785 | // Return the current volume "off" upward in the geometrical tree | |
786 | // ID and copy number | |
787 | if (off<0 || off>gGeoManager->GetLevel()) return 0; | |
788 | if (off==0) return CurrentVolID(copyNo); | |
789 | TGeoNode *node = gGeoManager->GetMother(off); | |
790 | if (!node) return 0; | |
791 | copyNo = node->GetNumber(); | |
2bc4c610 | 792 | if (fDebug) printf("CurrentVolOffId(off=%i,cpy=%i) = %i\n", off,copyNo,node->GetVolume()->GetNumber() ); |
8495a208 | 793 | return node->GetVolume()->GetNumber(); |
794 | } | |
795 | // FLUKA specific | |
796 | ||
797 | //_____________________________________________________________________________ | |
798 | const char* TFlukaMCGeometry::CurrentVolName() const | |
799 | { | |
800 | // | |
801 | // Returns the current volume name | |
802 | // | |
803 | if (gGeoManager->IsOutside()) return 0; | |
2bc4c610 | 804 | if (fDebug) printf("CurrentVolName : %s\n", gGeoManager->GetCurrentVolume()->GetName()); |
8495a208 | 805 | return gGeoManager->GetCurrentVolume()->GetName(); |
806 | } | |
807 | //_____________________________________________________________________________ | |
808 | const char* TFlukaMCGeometry::CurrentVolOffName(Int_t off) const | |
809 | { | |
810 | // | |
811 | // Return the current volume "off" upward in the geometrical tree | |
812 | // ID, name and copy number | |
813 | // if name=0 no name is returned | |
814 | // | |
815 | if (off<0 || off>gGeoManager->GetLevel()) return 0; | |
816 | if (off==0) return CurrentVolName(); | |
817 | TGeoNode *node = gGeoManager->GetMother(off); | |
818 | if (!node) return 0; | |
2bc4c610 | 819 | if (fDebug) printf("CurrentVolOffName(off=%i) : %s\n", off,node->GetVolume()->GetName()); |
8495a208 | 820 | return node->GetVolume()->GetName(); |
821 | } | |
822 | ||
823 | //_____________________________________________________________________________ | |
824 | void TFlukaMCGeometry::Gsatt(const char *name, const char *att, Int_t val) | |
825 | { | |
826 | // | |
827 | // NAME Volume name | |
828 | // IOPT Name of the attribute to be set | |
829 | // IVAL Value to which the attribute is to be set | |
830 | // see: TFluka::Gsatt | |
831 | char vname[5]; | |
832 | Vname(name,vname); | |
833 | char vatt[5]; | |
834 | Vname(att,vatt); | |
835 | gGeoManager->SetVolumeAttribute(vname, vatt, val); | |
836 | } | |
837 | ||
838 | //_____________________________________________________________________________ | |
839 | void TFlukaMCGeometry::Gdtom(Float_t *xd, Float_t *xm, Int_t iflag) | |
840 | { | |
841 | // | |
842 | // Computes coordinates XM (Master Reference System | |
843 | // knowing the coordinates XD (Detector Ref System) | |
844 | // The local reference system can be initialized by | |
845 | // - the tracking routines and GDTOM used in GUSTEP | |
846 | // - a call to GSCMED(NLEVEL,NAMES,NUMBER) | |
847 | // (inverse routine is GMTOD) | |
848 | // | |
849 | // If IFLAG=1 convert coordinates | |
850 | // IFLAG=2 convert direction cosinus | |
851 | // | |
852 | Double_t XM[3], XD[3]; | |
853 | Int_t i; | |
854 | for (i=0;i<3;i++) XD[i] = xd[i]; | |
855 | if (iflag == 1) gGeoManager->LocalToMaster(XD,XM); | |
856 | else gGeoManager->LocalToMasterVect(XD,XM); | |
857 | for (i=0;i<3;i++) xm[i]=XM[i]; | |
858 | } | |
859 | ||
860 | //_____________________________________________________________________________ | |
861 | void TFlukaMCGeometry::Gdtom(Double_t *xd, Double_t *xm, Int_t iflag) | |
862 | { | |
863 | if (iflag == 1) gGeoManager->LocalToMaster(xd,xm); | |
864 | else gGeoManager->LocalToMasterVect(xd,xm); | |
865 | } | |
866 | ||
867 | //_____________________________________________________________________________ | |
868 | void TFlukaMCGeometry::Gmtod(Float_t *xm, Float_t *xd, Int_t iflag) | |
869 | { | |
870 | // | |
871 | // Computes coordinates XD (in DRS) | |
872 | // from known coordinates XM in MRS | |
873 | // The local reference system can be initialized by | |
874 | // - the tracking routines and GMTOD used in GUSTEP | |
875 | // - a call to GMEDIA(XM,NUMED,CHECK) | |
876 | // - a call to GLVOLU(NLEVEL,NAMES,NUMBER,IER) | |
877 | // (inverse routine is GDTOM) | |
878 | // | |
879 | // If IFLAG=1 convert coordinates | |
880 | // IFLAG=2 convert direction cosinus | |
881 | // | |
882 | Double_t XM[3], XD[3]; | |
883 | Int_t i; | |
884 | for (i=0;i<3;i++) XM[i]=xm[i]; | |
885 | if (iflag == 1) gGeoManager->MasterToLocal(XM,XD); | |
886 | else gGeoManager->MasterToLocalVect(XM,XD); | |
887 | for (i=0;i<3;i++) xd[i] = XD[i]; | |
888 | } | |
889 | ||
890 | //_____________________________________________________________________________ | |
891 | void TFlukaMCGeometry::Gmtod(Double_t *xm, Double_t *xd, Int_t iflag) | |
892 | { | |
893 | if (iflag == 1) gGeoManager->MasterToLocal(xm,xd); | |
894 | else gGeoManager->MasterToLocalVect(xm,xd); | |
895 | } | |
896 | ||
897 | //_____________________________________________________________________________ | |
efde9b4d | 898 | void TFlukaMCGeometry::CreateFlukaMatFile(const char *fname) |
8495a208 | 899 | { |
900 | // ==== from FLUGG ==== | |
901 | // NAMES OF ELEMENTS AND COMPOUNDS: the names must be written in upper case, | |
902 | // according to the fluka standard. In addition,. they must be equal to the | |
903 | // names of the fluka materials - see fluka manual - in order that the | |
904 | // program load the right cross sections, and equal to the names included in | |
905 | // the .pemf. Otherwise the user must define the LOW-MAT CARDS, and make his | |
906 | // own .pemf, in order to get the right cross sections loaded in memory. | |
907 | ||
908 | Int_t zelem[128]; | |
909 | static char elNames[220] = { | |
910 | // 1 ============================= 5 ==================================== 10 ===================================== 15 === | |
911 | 'H','_','H','E','L','I','B','E','B','_','C','_','N','_','O','_','F','_','N','E','N','A','M','G','A','L','S','I','P','_', | |
912 | 'S','_','C','L','A','R','K','_','C','A','S','C','T','I','V','_','C','R','M','N','F','E','C','O','N','I','C','U','Z','N', | |
913 | 'G','A','G','E','A','S','S','E','B','R','K','R','R','B','S','R','Y','_','Z','R','N','B','M','O','T','C','R','U','R','H', | |
914 | 'P','D','A','G','C','D','I','N','S','N','S','B','T','E','I','_','X','E','C','S','B','A','L','A','C','E','P','R','N','D', | |
915 | 'P','M','S','M','E','U','G','D','T','B','D','Y','H','O','E','R','T','M','Y','B','L','U','H','F','T','A','W','_','R','E', | |
916 | 'O','S','I','R','P','T','A','U','H','G','T','L','P','B','B','I','P','O','A','T','R','N','F','R','R','A','A','C','T','H', | |
917 | 'P','A','U','_','N','P','P','U','A','M','C','M','B','K','C','F','E','S','F','M','M','D','N','O','L','R','R','F','D','B', | |
918 | 'S','G','B','H','H','S','M','T','D','S'}; | |
919 | memset(zelem, 0, 128*sizeof(Int_t)); | |
920 | TString sname; | |
921 | gGeoManager->Export("flgeom.root"); | |
922 | if (fname) sname = fname; | |
923 | else sname = "flukaMat.inp"; | |
924 | ofstream out; | |
925 | out.open(sname.Data(), ios::out); | |
926 | if (!out.good()) { | |
927 | Fatal("CreateFlukaMatFile", "could not open file %s for writing", sname.Data()); | |
928 | return; | |
929 | } | |
930 | PrintHeader(out, "MATERIALS AND COMPOUNDS"); | |
931 | PrintHeader(out, "MATERIALS"); | |
932 | TList *matlist = gGeoManager->GetListOfMaterials(); | |
933 | TIter next(matlist); | |
934 | Int_t nmater = matlist->GetSize(); | |
935 | Int_t nfmater = 0; | |
936 | TObjArray *listfluka = new TObjArray(nmater+50); | |
937 | TObjArray *listflukanames = new TObjArray(nmater+50); | |
938 | TGeoMaterial *mat, *matorig; | |
939 | TGeoMixture *mix = 0; | |
940 | TString matname; | |
941 | TObjString *objstr, *objstrother; | |
942 | Int_t i,j,k,idmat; | |
943 | Bool_t done; | |
944 | Int_t nelem, nidmat; | |
945 | Double_t amat,zmat,rhomat; | |
946 | Double_t zel, ael, wel, rho; | |
947 | char elname[8] = {' ',' ','_', 'E','L','E','M','\0'}; | |
948 | char digit[3]; | |
efde9b4d | 949 | Bool_t found = kFALSE; |
8495a208 | 950 | |
2bc4c610 | 951 | if (fDebug) printf("Creating materials and compounds\n"); |
8495a208 | 952 | for (i=0; i<nmater; i++) { |
953 | mat = (TGeoMaterial*)matlist->At(i); | |
efde9b4d | 954 | if (mat->GetZ()<1E-1) { |
955 | mat->SetIndex(2); // vacuum, built-in inside FLUKA | |
956 | continue; | |
957 | } | |
958 | // printf("material: %s index=%i: Z=%f A=%f rho=%f\n", mat->GetName(), mat->GetIndex(),mat->GetZ(),mat->GetA(),mat->GetDensity()); | |
8495a208 | 959 | matorig = gGeoManager->FindDuplicateMaterial(mat); |
960 | if (matorig) { | |
961 | idmat = matorig->GetIndex(); | |
962 | mat->SetIndex(idmat); | |
efde9b4d | 963 | // printf(" -> found a duplicate: %s with index %i\n", matorig->GetName(), idmat); |
8495a208 | 964 | matorig = 0; |
965 | } else { | |
efde9b4d | 966 | // printf(" Adding to temp list with index %i\n", nfmater+3); |
8495a208 | 967 | listfluka->Add(mat); |
968 | mat->SetIndex(nfmater+3); | |
969 | matorig = mat; | |
970 | objstr = new TObjString(mat->GetName()); | |
971 | listflukanames->Add(objstr); | |
972 | nfmater++; | |
973 | // look if name is existing | |
974 | nidmat = 0; | |
975 | matname = objstr->GetString(); | |
976 | ToFlukaString(matname); | |
977 | objstr->SetString(matname.Data()); | |
978 | done = kFALSE; | |
979 | while (!done) { | |
980 | if (nfmater == 1) break; | |
981 | for (j=0; j<nfmater-1; j++) { | |
982 | objstrother = (TObjString*)listflukanames->At(j); | |
983 | if (objstr->IsEqual(objstrother)) { | |
984 | // we have to change the name | |
985 | if (nidmat>98) { | |
986 | Error("CreateFlukaMatFile", "too many materials having same name"); | |
987 | return; | |
988 | } | |
989 | nidmat++; | |
990 | k = matname.Index(" "); | |
991 | if (k<0 || k>6) k=6; | |
992 | if (nidmat>9) { | |
993 | sprintf(digit, "%d", nidmat); | |
994 | } else { | |
995 | digit[0] = '0'; | |
996 | sprintf(&digit[1], "%d", nidmat); | |
997 | } | |
998 | matname.Insert(k,digit); | |
999 | matname.Remove(8); | |
1000 | objstr->SetString(matname.Data()); | |
1001 | break; | |
1002 | } | |
1003 | if (j == nfmater-2) { | |
1004 | done = kTRUE; | |
1005 | break; | |
1006 | } | |
1007 | } | |
1008 | } | |
efde9b4d | 1009 | // printf(" newmat name: %s\n", matname.Data()); |
8495a208 | 1010 | } |
1011 | // now we have unique materials with unique names in the lists | |
1012 | ||
8495a208 | 1013 | if (matorig && matorig->IsMixture()) { |
1014 | // create dummy materials for elements | |
1015 | rho = 0.999; | |
1016 | mix = (TGeoMixture*)matorig; | |
1017 | nelem = mix->GetNelements(); | |
efde9b4d | 1018 | // printf(" material is a MIXTURE with %i elements:\n", nelem); |
8495a208 | 1019 | for (j=0; j<nelem; j++) { |
efde9b4d | 1020 | found = kFALSE; |
8495a208 | 1021 | zel = (mix->GetZmixt())[j]; |
efde9b4d | 1022 | ael = (mix->GetAmixt())[j]; |
1023 | // printf(" Zelem[%i] = %g\n",j,zel); | |
8495a208 | 1024 | if ((zel-Int_t(zel))>0.01) { |
efde9b4d | 1025 | TGeoMaterial *mat1; |
1026 | for (Int_t imat=0; imat<nfmater; imat++) { | |
1027 | mat1 = (TGeoMaterial*)listfluka->At(imat); | |
1028 | if (TMath::Abs(mat1->GetZ()-zel)>1E-4) continue; | |
1029 | if (TMath::Abs(mat1->GetA()-ael)>1E-4) continue; | |
1030 | found = kTRUE; | |
1031 | break; | |
1032 | } | |
1033 | if (!found) Warning("CreateFlukaMatFile", "element with Z=%f\n", zel); | |
8495a208 | 1034 | } |
efde9b4d | 1035 | if (!zelem[Int_t(zel)] && !found) { |
8495a208 | 1036 | // write fluka element |
1037 | memcpy(elname, &elNames[2*Int_t(zel-1)], 2); | |
1038 | zelem[Int_t(zel)] = 1; | |
8495a208 | 1039 | mat = new TGeoMaterial(elname, ael, zel, rho); |
1040 | mat->SetIndex(nfmater+3); | |
efde9b4d | 1041 | // printf(" element not in list: new material %s at index=%i, Z=%g, A=%g, dummyrho=%g\n", |
1042 | // elname,nfmater+3,zel,ael,rho); | |
8495a208 | 1043 | listfluka->Add(mat); |
1044 | objstr = new TObjString(elname); | |
1045 | listflukanames->Add(objstr); | |
1046 | nfmater++; | |
1047 | } | |
1048 | } | |
1049 | } | |
1050 | } | |
1051 | // now dump materials in the file | |
efde9b4d | 1052 | // printf("DUMPING %i materials\n", nfmater); |
8495a208 | 1053 | for (i=0; i<nfmater; i++) { |
1054 | mat = (TGeoMaterial*)listfluka->At(i); | |
1055 | out << setw(10) << "MATERIAL "; | |
1056 | out.setf(static_cast<std::ios::fmtflags>(0),std::ios::floatfield); | |
b0853588 | 1057 | // matname = mat->GetName(); |
1058 | objstr = (TObjString*)listflukanames->At(i); | |
1059 | matname = objstr->GetString(); | |
8495a208 | 1060 | ToFlukaString(matname); |
1061 | zmat = mat->GetZ(); | |
efde9b4d | 1062 | if (zmat-Int_t(zmat)>0.01) { |
1063 | if (zmat-Int_t(zmat)>0.5) zmat = Int_t(zmat)+1.; | |
1064 | else zmat = Int_t(zmat); | |
1065 | } | |
8495a208 | 1066 | amat = mat->GetA(); |
1067 | rhomat = mat->GetDensity(); | |
1068 | // write material card | |
1069 | if (mat->IsMixture()) { | |
1070 | out << setw(10) << " "; | |
1071 | out << setw(10) << " "; | |
1072 | mix = (TGeoMixture*)mat; | |
1073 | } else { | |
1074 | out << setw(10) << setiosflags(ios::fixed) << setprecision(1) << zmat; | |
1075 | out << setw(10) << setprecision(3) << amat; | |
1076 | } | |
1077 | out.setf(static_cast<std::ios::fmtflags>(0),std::ios::floatfield); | |
1078 | out << setw(10) << setiosflags(ios::scientific) << setprecision(3) << rhomat; | |
1079 | out.setf(static_cast<std::ios::fmtflags>(0),std::ios::floatfield); | |
1080 | out << setw(10) << setiosflags(ios::fixed) << setprecision(1) << Double_t(i+3); | |
1081 | out << setw(10) << " "; | |
efde9b4d | 1082 | out << setw(10) << " "; |
1083 | out << setw(8) << matname.Data() << endl; | |
d11bf3ca | 1084 | // add LOW-MAT crd |
1085 | if (!mat->IsMixture()) { | |
1086 | out << setw(10) << "LOW-MAT "; | |
1087 | out.setf(static_cast<std::ios::fmtflags>(0),std::ios::floatfield); | |
1088 | out << setw(10) << setiosflags(ios::fixed) << setprecision(1) << Double_t(i+3); | |
1089 | out << setw(10) << " "; | |
1090 | out << setw(10) << " "; | |
1091 | out << setw(10) << " "; | |
1092 | out << setw(10) << " "; | |
1093 | out << setw(10) << " "; | |
1094 | out << setw(8) << matname.Data() << endl; | |
1095 | } | |
8495a208 | 1096 | } |
1097 | // write mixture header | |
1098 | PrintHeader(out, "COMPOUNDS"); | |
1099 | Int_t counttothree; | |
1100 | TGeoMaterial *element; | |
1101 | for (i=0; i<nfmater; i++) { | |
1102 | mat = (TGeoMaterial*)listfluka->At(i); | |
1103 | if (!mat->IsMixture()) continue; | |
1104 | mix = (TGeoMixture*)mat; | |
1105 | counttothree = 0; | |
1106 | out << setw(10) << "COMPOUND "; | |
1107 | nelem = mix->GetNelements(); | |
1108 | objstr = (TObjString*)listflukanames->At(i); | |
1109 | matname = objstr->GetString(); | |
efde9b4d | 1110 | // printf("MIXTURE %s with index %i having %i elements\n", matname.Data(), mat->GetIndex(),nelem); |
8495a208 | 1111 | for (j=0; j<nelem; j++) { |
1112 | // dump mixture cards | |
efde9b4d | 1113 | // printf(" #elem %i: Z=%g, A=%g, W=%g\n", j, (mix->GetZmixt())[j], |
1114 | // (mix->GetAmixt())[j],(mix->GetWmixt())[j]); | |
8495a208 | 1115 | wel = (mix->GetWmixt())[j]; |
1116 | zel = (mix->GetZmixt())[j]; | |
efde9b4d | 1117 | ael = (mix->GetAmixt())[j]; |
1118 | if (zel-Int_t(zel)>0.01) { | |
1119 | // loop the temporary list | |
1120 | element = 0; | |
1121 | TGeoMaterial *mat1; | |
1122 | for (Int_t imat=0; imat<i; imat++) { | |
1123 | mat1 = (TGeoMaterial*)listfluka->At(imat); | |
1124 | if (TMath::Abs(mat1->GetZ()-zel)>1E-4) continue; | |
1125 | if (TMath::Abs(mat1->GetA()-ael)>1E-4) continue; | |
1126 | element = mat1; | |
1127 | break; | |
1128 | } | |
1129 | } else { | |
1130 | memcpy(elname, &elNames[2*Int_t(zel-1)], 2); | |
1131 | element = (TGeoMaterial*)listfluka->FindObject(elname); | |
1132 | } | |
8495a208 | 1133 | if (!element) { |
1134 | Error("CreateFlukaMatFile", "Element Z=%g %s not found", zel, elname); | |
1135 | return; | |
1136 | } | |
1137 | idmat = element->GetIndex(); | |
efde9b4d | 1138 | // printf("element %s , index=%i\n", element->GetName(), idmat); |
8495a208 | 1139 | out.setf(static_cast<std::ios::fmtflags>(0),std::ios::floatfield); |
1140 | out << setw(10) << setiosflags(ios::fixed) << setprecision(6) << -wel; | |
1141 | out.setf(static_cast<std::ios::fmtflags>(0),std::ios::floatfield); | |
1142 | out << setw(10) << setiosflags(ios::fixed) << setprecision(1) << Double_t(idmat); | |
1143 | counttothree++; | |
1144 | if (counttothree == 3) { | |
1145 | out << matname.Data(); | |
1146 | out << endl; | |
1147 | if ( (j+1) != nelem) out << setw(10) << "COMPOUND "; | |
1148 | counttothree = 0; | |
1149 | } | |
1150 | } | |
1151 | //Unless we have 3, 6, 9... submaterials we need to put some empty | |
1152 | //space and the compound name | |
1153 | if (nelem%3) { | |
1154 | for (j=0; j<(3-(nelem%3)); j++) | |
1155 | out << setw(10) << " " << setw(10) << " "; | |
1156 | out << matname.Data(); | |
1157 | out << endl; | |
1158 | } | |
1159 | } | |
1160 | ||
1161 | // Now print the list of regions (volumes in TGeo) | |
1162 | Int_t nvols = gGeoManager->GetListOfUVolumes()->GetEntriesFast()-1; | |
1163 | TGeoVolume *vol; | |
1164 | /* | |
1165 | PrintHeader(out, "TGEO VOLUMES"); | |
1166 | for (i=1; i<=nvols; i++) { | |
1167 | vol = gGeoManager->GetVolume(i); | |
1168 | out.setf(std::ios::left, std::ios::adjustfield); | |
1169 | out << setw(10) << i; | |
1170 | out << setw(20) << vol->GetName() << endl; | |
1171 | } | |
1172 | */ | |
1173 | // Now print the material assignments | |
d11bf3ca | 1174 | Double_t flagfield = 0.; |
1175 | printf("#############################################################\n"); | |
1176 | if (fluka->IsFieldEnabled()) { | |
1177 | flagfield = 1.; | |
1178 | printf("Magnetic field enabled\n"); | |
1179 | } else printf("Magnetic field disabled\n"); | |
1180 | printf("#############################################################\n"); | |
1181 | ||
8495a208 | 1182 | PrintHeader(out, "TGEO MATERIAL ASSIGNMENTS"); |
1183 | for (i=1; i<=nvols; i++) { | |
1184 | vol = gGeoManager->GetVolume(i); | |
1185 | mat = vol->GetMedium()->GetMaterial(); | |
1186 | idmat = mat->GetIndex(); | |
8495a208 | 1187 | out << setw(10) << "ASSIGNMAT "; |
1188 | out.setf(static_cast<std::ios::fmtflags>(0),std::ios::floatfield); | |
1189 | out << setw(10) << setiosflags(ios::fixed) << Double_t(idmat); | |
1190 | out << setw(10) << setiosflags(ios::fixed) << Double_t(i); | |
1191 | out << setw(10) << "0.0"; | |
b0853588 | 1192 | out << setw(10) << "0.0"; |
8495a208 | 1193 | out << setw(10) << setiosflags(ios::fixed) << flagfield; |
b0853588 | 1194 | out << setw(10) << "0.0"; |
8495a208 | 1195 | out << endl; |
1196 | } | |
1197 | delete listfluka; | |
1198 | listflukanames->Delete(); | |
1199 | delete listflukanames; | |
1200 | out.close(); | |
efde9b4d | 1201 | fLastMaterial = nfmater+2; |
8495a208 | 1202 | } |
1203 | ||
1204 | //_____________________________________________________________________________ | |
1205 | void TFlukaMCGeometry::PrintHeader(ofstream &out, const char *text) const | |
1206 | { | |
1207 | // Print a FLUKA header. | |
1208 | out << "*\n" << "*\n" << "*\n"; | |
1209 | out << "********************* " << text << " *********************\n" | |
1210 | << "*\n"; | |
1211 | out << "*...+....1....+....2....+....3....+....4....+....5....+....6....+....7..." | |
1212 | << endl; | |
1213 | out << "*" << endl; | |
1214 | } | |
1215 | ||
1216 | //_____________________________________________________________________________ | |
1217 | Int_t TFlukaMCGeometry::RegionId() const | |
1218 | { | |
1219 | // Returns current region id <-> TGeo node id | |
1220 | if (gGeoManager->IsOutside()) return 0; | |
1221 | return gGeoManager->GetCurrentNode()->GetUniqueID(); | |
1222 | } | |
05265ca9 | 1223 | //_____________________________________________________________________________ |
1224 | void TFlukaMCGeometry::SetMreg(Int_t mreg) | |
1225 | { | |
1226 | // Update if needed next history; | |
d11bf3ca | 1227 | if (fluka->GetDummyBoundary()==2) { |
1228 | gGeoManager->CdNode(fNextLattice-1); | |
1229 | return; | |
1230 | } | |
05265ca9 | 1231 | Int_t curreg = (gGeoManager->IsOutside())?(mcgeom->NofVolumes()+1):gGeoManager->GetCurrentVolume()->GetNumber(); |
1232 | if (mreg==curreg) return; | |
1233 | if (mreg==fNextRegion) { | |
1234 | if (fNextLattice!=999999999) gGeoManager->CdNode(fNextLattice-1); | |
1235 | return; | |
d11bf3ca | 1236 | } else { |
1237 | if (mreg == fCurrentRegion) { | |
1238 | if (fCurrentLattice!=999999999) gGeoManager->CdNode(fCurrentLattice-1); | |
1239 | return; | |
1240 | } | |
1241 | } | |
2bc4c610 | 1242 | if (fDebug) printf("ERROR: mreg=%i neither current nor next region\n", mreg); |
05265ca9 | 1243 | } |
1244 | ||
d11bf3ca | 1245 | //_____________________________________________________________________________ |
1246 | void TFlukaMCGeometry::SetCurrentRegion(Int_t mreg, Int_t latt) | |
1247 | { | |
1248 | // Set index/history for next entered region | |
1249 | fCurrentRegion = mreg; | |
1250 | fCurrentLattice = latt; | |
1251 | } | |
1252 | ||
05265ca9 | 1253 | //_____________________________________________________________________________ |
1254 | void TFlukaMCGeometry::SetNextRegion(Int_t mreg, Int_t latt) | |
1255 | { | |
1256 | // Set index/history for next entered region | |
1257 | fNextRegion = mreg; | |
1258 | fNextLattice = latt; | |
1259 | } | |
8495a208 | 1260 | |
1261 | //_____________________________________________________________________________ | |
1262 | void TFlukaMCGeometry::ToFlukaString(TString &str) const | |
1263 | { | |
1264 | // ToFlukaString converts an string to something usefull in FLUKA: | |
1265 | // * Capital letters | |
1266 | // * Only 8 letters | |
1267 | // * Replace ' ' by '_' | |
1268 | if (str.Length()<8) { | |
1269 | str += " "; | |
1270 | } | |
1271 | str.Remove(8); | |
1272 | Int_t ilast; | |
1273 | for (ilast=7; ilast>0; ilast--) if (str(ilast)!=' ') break; | |
1274 | str.ToUpper(); | |
1275 | for (Int_t pos=0; pos<ilast; pos++) | |
1276 | if (str(pos)==' ') str.Replace(pos,1,"_",1); | |
1277 | return; | |
1278 | } | |
1279 | //______________________________________________________________________________ | |
1280 | void TFlukaMCGeometry::Vname(const char *name, char *vname) const | |
1281 | { | |
1282 | // | |
1283 | // convert name to upper case. Make vname at least 4 chars | |
1284 | // | |
1285 | Int_t l = strlen(name); | |
1286 | Int_t i; | |
1287 | l = l < 4 ? l : 4; | |
1288 | for (i=0;i<l;i++) vname[i] = toupper(name[i]); | |
1289 | for (i=l;i<4;i++) vname[i] = ' '; | |
1290 | vname[4] = 0; | |
1291 | } | |
1292 | ||
1293 | ||
1294 | // FLUKA GEOMETRY WRAPPERS - to replace FLUGG wrappers | |
1295 | ||
8495a208 | 1296 | //_____________________________________________________________________________ |
1297 | Int_t idnrwr(const Int_t & /*nreg*/, const Int_t & /*mlat*/) | |
1298 | { | |
1299 | // from FLUGG: | |
1300 | // Wrapper for setting DNEAR option on fluka side. Must return 0 | |
1301 | // if user doesn't want Fluka to use DNEAR to compute the | |
1302 | // step (the same effect is obtained with the GLOBAL (WHAT(3)=-1) | |
1303 | // card in fluka input), returns 1 if user wants Fluka always to | |
1304 | // use DNEAR (in this case, be sure that GEANT4 DNEAR is unique, | |
1305 | // coming from all directions!!!) | |
2bc4c610 | 1306 | if (mcgeom->IsDebugging()) printf("========== Dummy IDNRWR\n"); |
8495a208 | 1307 | return 0; |
1308 | } | |
1309 | ||
8495a208 | 1310 | //_____________________________________________________________________________ |
1311 | void g1wr(Double_t &pSx, Double_t &pSy, Double_t &pSz, | |
d11bf3ca | 1312 | Double_t *pV, Int_t &oldReg , const Int_t &oldLttc, Double_t &propStep, |
1313 | Int_t &nascFlag, Double_t &retStep, Int_t &newReg, | |
2573ac89 | 1314 | Double_t &saf, Int_t &newLttc, Int_t <tcFlag, |
8495a208 | 1315 | Double_t *sLt, Int_t *jrLt) |
d11bf3ca | 1316 | |
8495a208 | 1317 | { |
d11bf3ca | 1318 | // Initialize FLUKa point and direction; |
1319 | kNstep++; | |
1320 | /* | |
1321 | if (kNstep>0) { | |
1322 | mcgeom->SetDebugMode(kTRUE); | |
1323 | fluka->SetVerbosityLevel(3); | |
1324 | } | |
1325 | if (kNstep>6520) { | |
1326 | mcgeom->SetDebugMode(kFALSE); | |
1327 | fluka->SetVerbosityLevel(0); | |
1328 | } | |
1329 | if ((kNstep%10)==0) printf("step %i\n", kNstep); | |
1330 | */ | |
8495a208 | 1331 | |
2bc4c610 | 1332 | if (mcgeom->IsDebugging()) { |
1333 | printf("========== Inside G1WR\n"); | |
1334 | printf(" point/dir:(%14.9f, %14.9f, %14.9f, %g, %g, %g)\n", pSx,pSy,pSz,pV[0],pV[1],pV[2]); | |
d11bf3ca | 1335 | printf(" oldReg=%i oldLttc=%i pstep=%f\n",oldReg, oldLttc, propStep); |
2bc4c610 | 1336 | } |
8495a208 | 1337 | gGeoManager->SetCurrentPoint(pSx, pSy, pSz); |
1338 | gGeoManager->SetCurrentDirection(pV); | |
d11bf3ca | 1339 | mcgeom->SetCurrentRegion(oldReg, oldLttc); |
1340 | // Initialize default return values | |
1341 | lttcFlag = 0; | |
1342 | jrLt[lttcFlag] = oldLttc; | |
1343 | sLt[lttcFlag] = propStep; | |
1344 | jrLt[lttcFlag+1] = -1; | |
1345 | sLt[lttcFlag+1] = 0.; | |
1346 | newReg = oldReg; | |
1347 | newLttc = oldLttc; | |
1348 | // check if dummy boundary flag is set | |
1349 | Int_t curLttc, curReg; | |
1350 | if (fluka->IsDummyBoundary()) { | |
1351 | // printf("Dummy boundary intercepted. Point is: %f, %f, %f\n", pSx, pSy, pSz); | |
1352 | Bool_t crossedDummy = (oldLttc == TFlukaMCGeometry::kLttcVirtual)?kTRUE:kFALSE; | |
1353 | if (crossedDummy) { | |
1354 | // FLUKA crossed the dummy boundary - update new region/history | |
1355 | retStep = 0.; | |
1356 | saf = 0.; | |
1357 | mcgeom->GetNextRegion(newReg, newLttc); | |
1358 | mcgeom->SetMreg(newReg); | |
1359 | if (mcgeom->IsDebugging()) printf(" virtual newReg=%i newLttc=%i\n", newReg, newLttc); | |
1360 | sLt[lttcFlag] = 0.; // null step in current region | |
1361 | lttcFlag++; | |
1362 | jrLt[lttcFlag] = newLttc; | |
1363 | sLt[lttcFlag] = 0.; // null step in next region | |
1364 | jrLt[lttcFlag+1] = -1; | |
1365 | sLt[lttcFlag+1] = 0.; | |
1366 | fluka->SetDummyBoundary(0); | |
1367 | return; | |
1368 | } | |
1369 | } | |
1370 | ||
1371 | // Reset outside flag | |
05265ca9 | 1372 | if (gGeoManager->IsOutside()) { |
1373 | gGeoManager->SetOutside(kFALSE); | |
1374 | gGeoManager->CdTop(); | |
d11bf3ca | 1375 | } |
1376 | ||
1377 | // Reset dummy boundary flag | |
1378 | fluka->SetDummyBoundary(0); | |
1379 | ||
1380 | curLttc = gGeoManager->GetCurrentNodeId()+1; | |
1381 | curReg = gGeoManager->GetCurrentVolume()->GetNumber(); | |
2573ac89 | 1382 | if (oldLttc != curLttc) { |
d11bf3ca | 1383 | // FLUKA crossed the boundary : we trust that the given point is really there, |
1384 | // so we just update TGeo state | |
2573ac89 | 1385 | gGeoManager->CdNode(oldLttc-1); |
1386 | curLttc = gGeoManager->GetCurrentNodeId()+1; | |
d11bf3ca | 1387 | curReg = gGeoManager->GetCurrentVolume()->GetNumber(); |
1388 | if (mcgeom->IsDebugging()) printf(" re-initialized point: curReg=%i curLttc=%i\n", curReg, curLttc); | |
1389 | } | |
1390 | // Now the current TGeo state reflects the FLUKA state | |
1391 | if (mcgeom->IsDebugging()) printf(" current path: %s\n", gGeoManager->GetPath()); | |
1392 | Double_t extra = 1E-6; | |
1393 | Double_t tmpStep = propStep + extra; | |
1394 | gGeoManager->FindNextBoundary(-tmpStep); | |
1395 | Double_t snext = gGeoManager->GetStep(); | |
1396 | // !!!!! | |
1397 | if (snext<=0) { | |
1398 | // FLUKA is in the wrong region, notify it | |
1399 | if (mcgeom->IsDebugging()) printf("ERROR: snext=%f\n", snext); | |
1400 | // newReg = -3; | |
1401 | // return; | |
1402 | snext = extra; | |
1403 | } | |
1404 | saf = gGeoManager->GetSafeDistance(); | |
1405 | Bool_t cross = kFALSE; | |
1406 | Bool_t onBound = kFALSE; | |
1407 | if (snext<tmpStep) { | |
1408 | // We have some boundary in the way | |
1409 | Double_t dd = snext-propStep; | |
1410 | if (dd < 0) { | |
1411 | cross = kTRUE; | |
1412 | dd = -dd; | |
1413 | } | |
1414 | if (dd < 1E-8) onBound = kTRUE; | |
1415 | } | |
1416 | snext += 1.E-8; | |
1417 | if (mcgeom->IsDebugging()) { | |
1418 | if (!cross) printf(" physical step approved: %f\n", propStep); | |
1419 | else printf(" boundary crossing at: %f\n", snext); | |
1420 | if (onBound) printf(" step on boundary limit ! NASC=%i\n", nascFlag); | |
1421 | } | |
1422 | if (!cross) { | |
1423 | // Next boundary further than proposed step, which is approved | |
1424 | retStep = propStep; | |
1425 | sLt[lttcFlag] = propStep; | |
1426 | return; | |
1427 | } | |
1428 | // The next boundary is closer. We try to cross it. | |
2573ac89 | 1429 | Double_t *point = gGeoManager->GetCurrentPoint(); |
1430 | Double_t *dir = gGeoManager->GetCurrentDirection(); | |
d11bf3ca | 1431 | Double_t pt[3]; |
1432 | memcpy(pt, point, 3*sizeof(Double_t)); | |
1433 | ||
2573ac89 | 1434 | Int_t i; |
d11bf3ca | 1435 | for (i=0;i<3;i++) point[i] += snext*dir[i]; |
1436 | gGeoManager->FindNode(); | |
1437 | newLttc = (gGeoManager->IsOutside())?(TFlukaMCGeometry::kLttcOutside):gGeoManager->GetCurrentNodeId()+1; | |
1438 | if (newLttc == oldLttc) { | |
1439 | // brute force ... | |
1440 | // Just try a fast extra small step | |
1441 | snext += 1E-6; | |
1442 | for (i=0;i<3;i++) point[i] = pt[i]+snext*dir[i]; | |
1443 | gGeoManager->FindNode(); | |
1444 | newLttc = (gGeoManager->IsOutside())?(TFlukaMCGeometry::kLttcOutside):gGeoManager->GetCurrentNodeId()+1; | |
1445 | if (newLttc == oldLttc) { | |
1446 | // check if state changes at the end of the proposed step | |
1447 | for (i=0;i<3;i++) point[i] = pt[i]+propStep*dir[i]; | |
2573ac89 | 1448 | gGeoManager->FindNode(); |
d11bf3ca | 1449 | newLttc = (gGeoManager->IsOutside())?(TFlukaMCGeometry::kLttcOutside):gGeoManager->GetCurrentNodeId()+1; |
1450 | if (newLttc==oldLttc) { | |
1451 | // approve step | |
1452 | retStep = propStep; | |
1453 | sLt[lttcFlag] = propStep; | |
1454 | return; | |
05265ca9 | 1455 | } |
d11bf3ca | 1456 | // snext is underestimated - we will create a virtual one to overcome the error |
1457 | // printf("some boundary in the way...\n"); | |
1458 | } | |
1459 | } | |
1460 | gGeoManager->SetCurrentPoint(pt); | |
1461 | // newLttc = (gGeoManager->IsOutside())?(TFlukaMCGeometry::kLttcOutside):gGeoManager->GetCurrentNodeId()+1; | |
1462 | newReg = (gGeoManager->IsOutside())?(mcgeom->NofVolumes()+1):gGeoManager->GetCurrentVolume()->GetNumber(); | |
1463 | if (mcgeom->IsDebugging()) printf(" newReg=%i newLttc=%i\n", newReg, newLttc); | |
1464 | ||
1465 | // We really crossed the boundary, but is it the same region ? | |
05265ca9 | 1466 | mcgeom->SetNextRegion(newReg, newLttc); |
d11bf3ca | 1467 | if (newReg == oldReg) { |
1468 | // Virtual boundary between replicants | |
1469 | if (mcgeom->IsDebugging()) printf(" DUMMY boundary\n"); | |
1470 | newReg = 1; // cheat FLUKA telling it it crossed the TOP region | |
1471 | newLttc = TFlukaMCGeometry::kLttcVirtual; | |
1472 | // mark that next boundary is virtual | |
1473 | fluka->SetDummyBoundary(1); | |
1474 | } | |
1475 | retStep = snext; | |
1476 | sLt[lttcFlag] = snext; | |
1477 | lttcFlag++; | |
1478 | jrLt[lttcFlag] = newLttc; | |
1479 | sLt[lttcFlag] = snext; | |
1480 | jrLt[lttcFlag+1] = -1; | |
1481 | sLt[lttcFlag+1] = 0.; | |
1482 | ||
2573ac89 | 1483 | if (newLttc!=oldLttc) { |
1484 | if (gGeoManager->IsOutside()) { | |
1485 | gGeoManager->SetOutside(kFALSE); | |
1486 | gGeoManager->CdTop(); | |
d11bf3ca | 1487 | } |
1488 | gGeoManager->CdTop(); | |
1489 | if (!gGeoManager->GetCurrentMatrix()->IsIdentity()) printf("ERROR at step %i\n", kNstep); | |
2573ac89 | 1490 | gGeoManager->CdNode(oldLttc-1); |
8495a208 | 1491 | } |
d11bf3ca | 1492 | if (mcgeom->IsDebugging()) { |
1493 | printf("=> snext=%g safe=%g\n", snext, saf); | |
1494 | for (Int_t i=0; i<lttcFlag+1; i++) printf(" jrLt[%i]=%i sLt[%i]=%g\n", i,jrLt[i],i,sLt[i]); | |
1495 | } | |
2bc4c610 | 1496 | if (mcgeom->IsDebugging()) printf("<= G1WR (in: %s)\n", gGeoManager->GetPath()); |
8495a208 | 1497 | } |
1498 | ||
efde9b4d | 1499 | //_____________________________________________________________________________ |
1500 | void g1rtwr() | |
1501 | { | |
2bc4c610 | 1502 | if (mcgeom->IsDebugging()) printf("========== Dummy G1RTWR\n"); |
efde9b4d | 1503 | } |
1504 | ||
1505 | //_____________________________________________________________________________ | |
1506 | void conhwr(Int_t & /*intHist*/, Int_t * /*incrCount*/) | |
1507 | { | |
2bc4c610 | 1508 | if (mcgeom->IsDebugging()) printf("========== Dummy CONHWR\n"); |
efde9b4d | 1509 | } |
1510 | ||
1511 | //_____________________________________________________________________________ | |
2573ac89 | 1512 | void inihwr(Int_t &intHist) |
efde9b4d | 1513 | { |
2bc4c610 | 1514 | if (mcgeom->IsDebugging()) printf("========== Inside INIHWR -> reinitializing history: %i\n", intHist); |
2573ac89 | 1515 | if (gGeoManager->IsOutside()) gGeoManager->CdTop(); |
1516 | if (intHist<=0) { | |
1517 | // printf("=== wrong history number\n"); | |
1518 | return; | |
1519 | } | |
1520 | if (intHist==0) gGeoManager->CdTop(); | |
1521 | else gGeoManager->CdNode(intHist-1); | |
2bc4c610 | 1522 | if (mcgeom->IsDebugging()) { |
1523 | printf(" --- current path: %s\n", gGeoManager->GetPath()); | |
1524 | printf("<= INIHWR\n"); | |
1525 | } | |
efde9b4d | 1526 | } |
1527 | ||
1528 | //_____________________________________________________________________________ | |
1529 | void jomiwr(const Int_t & /*nge*/, const Int_t & /*lin*/, const Int_t & /*lou*/, | |
1530 | Int_t &flukaReg) | |
1531 | { | |
1532 | // Geometry initialization wrapper called by FLUKAM. Provides to FLUKA the | |
1533 | // number of regions (volumes in TGeo) | |
1534 | // build application geometry | |
2bc4c610 | 1535 | if (mcgeom->IsDebugging()) printf("========== Inside JOMIWR\n"); |
2573ac89 | 1536 | flukaReg = gGeoManager->GetListOfUVolumes()->GetEntriesFast(); |
2bc4c610 | 1537 | if (mcgeom->IsDebugging()) printf("<= JOMIWR: last region=%i\n", flukaReg); |
efde9b4d | 1538 | } |
1539 | ||
1540 | //_____________________________________________________________________________ | |
1541 | void lkdbwr(Double_t &pSx, Double_t &pSy, Double_t &pSz, | |
1542 | Double_t * /*pV*/, const Int_t &oldReg, const Int_t &oldLttc, | |
1543 | Int_t &newReg, Int_t &flagErr, Int_t &newLttc) | |
1544 | { | |
2bc4c610 | 1545 | if (mcgeom->IsDebugging()) { |
1546 | printf("========== Inside LKDBWR (%f, %f, %f)\n",pSx, pSy, pSz); | |
1547 | // printf(" in: pV=(%f, %f, %f)\n", pV[0], pV[1], pV[2]); | |
1548 | printf(" in: oldReg=%i oldLttc=%i\n", oldReg, oldLttc); | |
1549 | } | |
efde9b4d | 1550 | TGeoNode *node = gGeoManager->FindNode(pSx, pSy, pSz); |
1551 | if (gGeoManager->IsOutside()) { | |
1552 | newReg = mcgeom->NofVolumes()+1; | |
05265ca9 | 1553 | // newLttc = gGeoManager->GetCurrentNodeId(); |
1554 | newLttc = 999999999; | |
2bc4c610 | 1555 | if (mcgeom->IsDebugging()) { |
1556 | printf("OUTSIDE\n"); | |
1557 | printf(" out: newReg=%i newLttc=%i\n", newReg, newLttc); | |
1558 | printf("<= LKMGWR\n"); | |
1559 | } | |
05265ca9 | 1560 | flagErr = newReg; |
1561 | return; | |
efde9b4d | 1562 | } |
1563 | newReg = node->GetVolume()->GetNumber(); | |
1564 | newLttc = gGeoManager->GetCurrentNodeId()+1; | |
d11bf3ca | 1565 | mcgeom->SetNextRegion(newReg, newLttc); |
efde9b4d | 1566 | flagErr = newReg; |
2bc4c610 | 1567 | if (mcgeom->IsDebugging()) { |
1568 | printf(" out: newReg=%i newLttc=%i\n", newReg, newLttc); | |
1569 | printf("<= LKDBWR\n"); | |
1570 | } | |
efde9b4d | 1571 | } |
1572 | ||
1573 | //_____________________________________________________________________________ | |
1574 | void lkfxwr(Double_t &pSx, Double_t &pSy, Double_t &pSz, | |
1575 | Double_t * /*pV*/, const Int_t &oldReg, const Int_t &oldLttc, | |
1576 | Int_t &newReg, Int_t &flagErr, Int_t &newLttc) | |
1577 | { | |
2bc4c610 | 1578 | if (mcgeom->IsDebugging()) { |
1579 | printf("========== Inside LKFXWR (%f, %f, %f)\n",pSx, pSy, pSz); | |
1580 | // printf(" in: pV=(%f, %f, %f)\n", pV[0], pV[1], pV[2]); | |
1581 | printf(" in: oldReg=%i oldLttc=%i\n", oldReg, oldLttc); | |
1582 | } | |
efde9b4d | 1583 | TGeoNode *node = gGeoManager->FindNode(pSx, pSy, pSz); |
1584 | if (gGeoManager->IsOutside()) { | |
1585 | newReg = mcgeom->NofVolumes()+1; | |
05265ca9 | 1586 | // newLttc = gGeoManager->GetCurrentNodeId(); |
1587 | newLttc = 999999999; | |
2bc4c610 | 1588 | if (mcgeom->IsDebugging()) { |
1589 | printf("OUTSIDE\n"); | |
1590 | printf(" out: newReg=%i newLttc=%i\n", newReg, newLttc); | |
1591 | printf("<= LKMGWR\n"); | |
1592 | } | |
05265ca9 | 1593 | flagErr = newReg; |
1594 | return; | |
efde9b4d | 1595 | } |
1596 | newReg = node->GetVolume()->GetNumber(); | |
1597 | newLttc = gGeoManager->GetCurrentNodeId()+1; | |
d11bf3ca | 1598 | mcgeom->SetNextRegion(newReg, newLttc); |
efde9b4d | 1599 | flagErr = newReg; |
2bc4c610 | 1600 | if (mcgeom->IsDebugging()) { |
1601 | printf(" out: newReg=%i newLttc=%i\n", newReg, newLttc); | |
1602 | printf("<= LKFXWR\n"); | |
1603 | } | |
efde9b4d | 1604 | } |
1605 | ||
1606 | //_____________________________________________________________________________ | |
1607 | void lkmgwr(Double_t &pSx, Double_t &pSy, Double_t &pSz, | |
1608 | Double_t * /*pV*/, const Int_t &oldReg, const Int_t &oldLttc, | |
1609 | Int_t &flagErr, Int_t &newReg, Int_t &newLttc) | |
1610 | { | |
2bc4c610 | 1611 | if (mcgeom->IsDebugging()) { |
1612 | printf("========== Inside LKMGWR (%f, %f, %f)\n",pSx, pSy, pSz); | |
1613 | // printf(" in: pV=(%f, %f, %f)\n", pV[0], pV[1], pV[2]); | |
1614 | printf(" in: oldReg=%i oldLttc=%i\n", oldReg, oldLttc); | |
1615 | } | |
efde9b4d | 1616 | TGeoNode *node = gGeoManager->FindNode(pSx, pSy, pSz); |
1617 | if (gGeoManager->IsOutside()) { | |
1618 | newReg = mcgeom->NofVolumes()+1; | |
05265ca9 | 1619 | // newLttc = gGeoManager->GetCurrentNodeId(); |
1620 | newLttc = 999999999; | |
2bc4c610 | 1621 | if (mcgeom->IsDebugging()) { |
1622 | printf("OUTSIDE\n"); | |
1623 | printf(" out: newReg=%i newLttc=%i\n", newReg, newLttc); | |
1624 | printf("<= LKMGWR\n"); | |
1625 | } | |
05265ca9 | 1626 | flagErr = newReg; |
1627 | return; | |
efde9b4d | 1628 | } |
1629 | newReg = node->GetVolume()->GetNumber(); | |
1630 | newLttc = gGeoManager->GetCurrentNodeId()+1; | |
d11bf3ca | 1631 | mcgeom->SetNextRegion(newReg, newLttc); |
efde9b4d | 1632 | flagErr = newReg; |
2bc4c610 | 1633 | if (mcgeom->IsDebugging()) { |
1634 | printf(" out: newReg=%i newLttc=%i\n", newReg, newLttc); | |
1635 | printf("<= LKMGWR\n"); | |
1636 | } | |
efde9b4d | 1637 | } |
1638 | ||
1639 | //_____________________________________________________________________________ | |
1640 | void lkwr(Double_t &pSx, Double_t &pSy, Double_t &pSz, | |
1641 | Double_t * /*pV*/, const Int_t &oldReg, const Int_t &oldLttc, | |
1642 | Int_t &newReg, Int_t &flagErr, Int_t &newLttc) | |
1643 | { | |
2bc4c610 | 1644 | if (mcgeom->IsDebugging()) { |
1645 | printf("========== Inside LKWR (%f, %f, %f)\n",pSx, pSy, pSz); | |
1646 | // printf(" in: pV=(%f, %f, %f)\n", pV[0], pV[1], pV[2]); | |
1647 | printf(" in: oldReg=%i oldLttc=%i\n", oldReg, oldLttc); | |
1648 | } | |
efde9b4d | 1649 | TGeoNode *node = gGeoManager->FindNode(pSx, pSy, pSz); |
1650 | if (gGeoManager->IsOutside()) { | |
1651 | newReg = mcgeom->NofVolumes()+1; | |
05265ca9 | 1652 | // newLttc = gGeoManager->GetCurrentNodeId(); |
1653 | newLttc = 999999999; | |
2bc4c610 | 1654 | if (mcgeom->IsDebugging()) { |
1655 | printf("OUTSIDE\n"); | |
1656 | printf(" out: newReg=%i newLttc=%i\n", newReg, newLttc); | |
1657 | printf("<= LKMGWR\n"); | |
1658 | } | |
05265ca9 | 1659 | flagErr = newReg; |
1660 | return; | |
efde9b4d | 1661 | } |
1662 | newReg = node->GetVolume()->GetNumber(); | |
d11bf3ca | 1663 | newLttc = gGeoManager->GetCurrentNodeId()+1; |
1664 | mcgeom->SetNextRegion(newReg, newLttc); | |
efde9b4d | 1665 | flagErr = newReg; |
2bc4c610 | 1666 | if (mcgeom->IsDebugging()) { |
1667 | printf(" out: newReg=%i newLttc=%i in %s\n", newReg, newLttc, gGeoManager->GetPath()); | |
1668 | printf("<= LKWR\n"); | |
1669 | } | |
efde9b4d | 1670 | } |
1671 | ||
1672 | //_____________________________________________________________________________ | |
2573ac89 | 1673 | void nrmlwr(Double_t &pSx, Double_t &pSy, Double_t &pSz, |
1674 | Double_t &pVx, Double_t &pVy, Double_t &pVz, | |
1675 | Double_t *norml, const Int_t &oldReg, | |
1676 | const Int_t &newReg, Int_t &flagErr) | |
efde9b4d | 1677 | { |
2bc4c610 | 1678 | if (mcgeom->IsDebugging()) { |
1679 | printf("========== Inside NRMLWR (%g, %g, %g, %g, %g, %g)\n", pSx,pSy,pSz,pVx,pVy,pVz); | |
1680 | printf(" oldReg=%i, newReg=%i\n", oldReg,newReg); | |
1681 | } | |
d11bf3ca | 1682 | // Int_t curreg = (gGeoManager->IsOutside())?(mcgeom->NofVolumes()+1):gGeoManager->GetCurrentVolume()->GetNumber(); |
1683 | // Int_t curLttc = gGeoManager->GetCurrentNodeId()+1; | |
1684 | // if (mcgeom->IsDebugging()) printf(" curReg=%i, curLttc=%i in: %s\n", curreg, curLttc, gGeoManager->GetPath()); | |
1685 | // Bool_t regsame = (curreg==oldReg)?kTRUE:kFALSE; | |
2573ac89 | 1686 | gGeoManager->SetCurrentPoint(pSx, pSy, pSz); |
1687 | gGeoManager->SetCurrentDirection(pVx,pVy,pVz); | |
d11bf3ca | 1688 | /* |
2573ac89 | 1689 | if (!regsame) { |
2bc4c610 | 1690 | if (mcgeom->IsDebugging()) printf(" REGIONS DOEN NOT MATCH\n"); |
2573ac89 | 1691 | gGeoManager->FindNode(); |
1692 | curreg = (gGeoManager->IsOutside())?(mcgeom->NofVolumes()+1):gGeoManager->GetCurrentVolume()->GetNumber(); | |
1693 | curLttc = gGeoManager->GetCurrentNodeId()+1; | |
2bc4c610 | 1694 | if (mcgeom->IsDebugging()) printf(" re-initialized point: curReg=%i curLttc=%i curPath=%s\n", curreg, curLttc, gGeoManager->GetPath()); |
2573ac89 | 1695 | } |
d11bf3ca | 1696 | */ |
2573ac89 | 1697 | Double_t *dnorm = gGeoManager->FindNormalFast(); |
1698 | flagErr = 0; | |
1699 | if (!dnorm) { | |
1700 | printf(" ERROR: Cannot compute fast normal\n"); | |
1701 | flagErr = 1; | |
1702 | norml[0] = -pVx; | |
1703 | norml[1] = -pVy; | |
1704 | norml[2] = -pVz; | |
1705 | } | |
1706 | norml[0] = -dnorm[0]; | |
1707 | norml[1] = -dnorm[1]; | |
1708 | norml[2] = -dnorm[2]; | |
2bc4c610 | 1709 | if (mcgeom->IsDebugging()) printf(" normal to boundary: (%g, %g, %g)\n", norml[0], norml[1], norml[2]); |
d11bf3ca | 1710 | // curreg = (gGeoManager->IsOutside())?(mcgeom->NofVolumes()+1):gGeoManager->GetCurrentVolume()->GetNumber(); |
1711 | // curLttc = gGeoManager->GetCurrentNodeId()+1; | |
2bc4c610 | 1712 | if (mcgeom->IsDebugging()) { |
d11bf3ca | 1713 | // printf(" final location: curReg=%i, curLttc=%i in %s\n", curreg,curLttc,gGeoManager->GetPath()); |
2bc4c610 | 1714 | printf("<= NRMLWR\n"); |
1715 | } | |
efde9b4d | 1716 | } |
1717 | ||
1718 | //_____________________________________________________________________________ | |
1719 | void rgrpwr(const Int_t & /*flukaReg*/, const Int_t & /*ptrLttc*/, Int_t & /*g4Reg*/, | |
1720 | Int_t * /*indMother*/, Int_t * /*repMother*/, Int_t & /*depthFluka*/) | |
1721 | { | |
2bc4c610 | 1722 | if (mcgeom->IsDebugging()) printf("=> Dummy RGRPWR\n"); |
efde9b4d | 1723 | } |
1724 | ||
1725 | //_____________________________________________________________________________ | |
1726 | Int_t isvhwr(const Int_t &check, const Int_t & intHist) | |
1727 | { | |
1728 | // from FLUGG: | |
1729 | // Wrapper for saving current navigation history (fCheck=default) | |
1730 | // and returning its pointer. If fCheck=-1 copy of history pointed | |
1731 | // by intHist is made in NavHistWithCount object, and its pointer | |
1732 | // is returned. fCheck=1 and fCheck=2 cases are only in debugging | |
1733 | // version: an array is created by means of FGeometryInit functions | |
1734 | // (but could be a static int * ptrArray = new int[10000] with | |
1735 | // file scope as well) that stores a flag for deleted/undeleted | |
1736 | // histories and at the end of event is checked to verify that | |
1737 | // all saved history objects have been deleted. | |
1738 | ||
1739 | // For TGeo, just return the current node ID. No copy need to be made. | |
1740 | ||
2bc4c610 | 1741 | if (mcgeom->IsDebugging()) printf("=> Inside ISVHWR\n"); |
efde9b4d | 1742 | if (check<0) return intHist; |
2573ac89 | 1743 | Int_t histInt = gGeoManager->GetCurrentNodeId()+1; |
2bc4c610 | 1744 | if (mcgeom->IsDebugging()) printf("<= ISVHWR: history is: %i in: %s\n", histInt, gGeoManager->GetPath()); |
efde9b4d | 1745 | return histInt; |
1746 | } | |
1747 | ||
1748 | ||
8495a208 | 1749 | |
1750 |