]>
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" |
c1a0a2f0 | 27 | #include "TCallf77.h" |
89125dc7 | 28 | #include "TFluka.h" |
8495a208 | 29 | #include "TFlukaMCGeometry.h" |
6d184c54 | 30 | #include "TFlukaConfigOption.h" |
8495a208 | 31 | #include "TGeoManager.h" |
32 | #include "TGeoVolume.h" | |
c1a0a2f0 | 33 | #include "TObjString.h" |
81f1d030 | 34 | #include "Fsourcm.h" |
3860051a | 35 | |
8495a208 | 36 | #ifndef WIN32 |
37 | # define idnrwr idnrwr_ | |
38 | # define g1wr g1wr_ | |
39 | # define g1rtwr g1rtwr_ | |
40 | # define conhwr conhwr_ | |
41 | # define inihwr inihwr_ | |
42 | # define jomiwr jomiwr_ | |
43 | # define lkdbwr lkdbwr_ | |
44 | # define lkfxwr lkfxwr_ | |
45 | # define lkmgwr lkmgwr_ | |
46 | # define lkwr lkwr_ | |
47 | # define magfld magfld_ | |
48 | # define nrmlwr nrmlwr_ | |
49 | # define rgrpwr rgrpwr_ | |
50 | # define isvhwr isvhwr_ | |
51 | ||
52 | #else | |
53 | ||
54 | # define idnrwr IDNRWR | |
55 | # define g1wr G1WR | |
56 | # define g1rtwr G1RTWR | |
57 | # define conhwr CONHWR | |
58 | # define inihwr INIHWR | |
59 | # define jomiwr JOMIWR | |
60 | # define lkdbwr LKDBWR | |
61 | # define lkfxwr LKFXWR | |
62 | # define lkmgwr LKMGWR | |
63 | # define lkwr LKWR | |
64 | # define magfld MAGFLD | |
65 | # define nrmlwr NRMLWR | |
66 | # define rgrpwr RGRPWR | |
67 | # define isvhwr ISVHWR | |
68 | ||
69 | #endif | |
70 | ||
71 | //____________________________________________________________________________ | |
72 | extern "C" | |
73 | { | |
74 | // | |
75 | // Prototypes for FLUKA navigation methods | |
76 | // | |
77 | Int_t type_of_call idnrwr(const Int_t & /*nreg*/, const Int_t & /*mlat*/); | |
78 | void type_of_call g1wr(Double_t & /*pSx*/, Double_t & /*pSy*/, Double_t & /*pSz*/, | |
79 | Double_t * /*pV*/, Int_t & /*oldReg*/ , const Int_t & /*oldLttc*/, Double_t & /*propStep*/, | |
80 | Int_t & /*nascFlag*/, Double_t & /*retStep*/, Int_t & /*newReg*/, | |
81 | Double_t & /*saf*/, Int_t & /*newLttc*/, Int_t & /*LttcFlag*/, | |
82 | Double_t *s /*Lt*/, Int_t * /*jrLt*/); | |
83 | ||
84 | void type_of_call g1rtwr(); | |
85 | void type_of_call conhwr(Int_t & /*intHist*/, Int_t * /*incrCount*/); | |
86 | void type_of_call inihwr(Int_t & /*intHist*/); | |
87 | void type_of_call jomiwr(const Int_t & /*nge*/, const Int_t & /*lin*/, const Int_t & /*lou*/, | |
88 | Int_t & /*flukaReg*/); | |
89 | void type_of_call lkdbwr(Double_t & /*pSx*/, Double_t & /*pSy*/, Double_t & /*pSz*/, | |
90 | Double_t * /*pV*/, const Int_t & /*oldReg*/, const Int_t & /*oldLttc*/, | |
91 | Int_t & /*newReg*/, Int_t & /*flagErr*/, Int_t & /*newLttc*/); | |
92 | void type_of_call lkfxwr(Double_t & /*pSx*/, Double_t & /*pSy*/, Double_t & /*pSz*/, | |
93 | Double_t * /*pV*/, const Int_t & /*oldReg*/, const Int_t & /*oldLttc*/, | |
94 | Int_t & /*newReg*/, Int_t & /*flagErr*/, Int_t & /*newLttc*/); | |
95 | void type_of_call lkmgwr(Double_t & /*pSx*/, Double_t & /*pSy*/, Double_t & /*pSz*/, | |
96 | Double_t * /*pV*/, const Int_t & /*oldReg*/, const Int_t & /*oldLttc*/, | |
97 | Int_t & /*flagErr*/, Int_t & /*newReg*/, Int_t & /*newLttc*/); | |
98 | void type_of_call lkwr(Double_t & /*pSx*/, Double_t & /*pSy*/, Double_t & /*pSz*/, | |
99 | Double_t * /*pV*/, const Int_t & /*oldReg*/, const Int_t & /*oldLttc*/, | |
100 | Int_t & /*newReg*/, Int_t & /*flagErr*/, Int_t & /*newLttc*/); | |
efde9b4d | 101 | // void type_of_call magfld(const Double_t & /*pX*/, const Double_t & /*pY*/, const Double_t & /*pZ*/, |
102 | // Double_t & /*cosBx*/, Double_t & /*cosBy*/, Double_t & /*cosBz*/, | |
103 | // Double_t & /*Bmag*/, Int_t & /*reg*/, Int_t & /*idiscflag*/); | |
8495a208 | 104 | void type_of_call nrmlwr(Double_t & /*pSx*/, Double_t & /*pSy*/, Double_t & /*pSz*/, |
105 | Double_t & /*pVx*/, Double_t & /*pVy*/, Double_t & /*pVz*/, | |
106 | Double_t * /*norml*/, const Int_t & /*oldReg*/, | |
107 | const Int_t & /*newReg*/, Int_t & /*flagErr*/); | |
108 | void type_of_call rgrpwr(const Int_t & /*flukaReg*/, const Int_t & /*ptrLttc*/, Int_t & /*g4Reg*/, | |
109 | Int_t * /*indMother*/, Int_t * /*repMother*/, Int_t & /*depthFluka*/); | |
110 | Int_t type_of_call isvhwr(const Int_t & /*fCheck*/, const Int_t & /*intHist*/); | |
8be76537 | 111 | } |
8495a208 | 112 | |
113 | // TFluka global pointer | |
89125dc7 | 114 | TFluka *gFluka = 0; |
115 | TFlukaMCGeometry *gMCGeom = 0; | |
6f5667d1 | 116 | Int_t gNstep = 0; |
8495a208 | 117 | |
118 | ClassImp(TFlukaMCGeometry) | |
119 | ||
81f1d030 | 120 | TFlukaMCGeometry* TFlukaMCGeometry::fgInstance= NULL; |
8495a208 | 121 | |
122 | //_____________________________________________________________________________ | |
123 | TFlukaMCGeometry::TFlukaMCGeometry(const char *name, const char *title) | |
89125dc7 | 124 | : TNamed(name, title) |
8495a208 | 125 | { |
126 | // | |
127 | // Standard constructor | |
128 | // | |
2bc4c610 | 129 | fDebug = kFALSE; |
efde9b4d | 130 | fLastMaterial = 0; |
81f1d030 | 131 | fCurrentRegion = 0; |
132 | fCurrentLattice = 0; | |
b1536e91 | 133 | fNextRegion = 0; |
134 | fNextLattice = 0; | |
135 | fRegionList = 0; | |
81f1d030 | 136 | fIndmat = 0; |
89125dc7 | 137 | gFluka = (TFluka*)gMC; |
138 | gMCGeom = this; | |
6f5667d1 | 139 | gNstep = 0; |
89125dc7 | 140 | fMatList = new TObjArray(256); |
141 | fMatNames = new TObjArray(256); | |
8495a208 | 142 | } |
143 | ||
144 | //_____________________________________________________________________________ | |
145 | TFlukaMCGeometry::TFlukaMCGeometry() | |
8495a208 | 146 | { |
147 | // | |
148 | // Default constructor | |
149 | // | |
2bc4c610 | 150 | fDebug = kFALSE; |
efde9b4d | 151 | fLastMaterial = 0; |
81f1d030 | 152 | fCurrentRegion = 0; |
153 | fCurrentLattice = 0; | |
b1536e91 | 154 | fNextRegion = 0; |
155 | fNextLattice = 0; | |
156 | fRegionList = 0; | |
81f1d030 | 157 | fIndmat = 0; |
89125dc7 | 158 | gFluka = (TFluka*)gMC; |
159 | gMCGeom = this; | |
6f5667d1 | 160 | gNstep = 0; |
89125dc7 | 161 | fMatList = 0; |
162 | fMatNames = 0; | |
8495a208 | 163 | } |
164 | ||
165 | //_____________________________________________________________________________ | |
166 | TFlukaMCGeometry::~TFlukaMCGeometry() | |
167 | { | |
168 | // | |
169 | // Destructor | |
170 | // | |
171 | fgInstance=0; | |
b1536e91 | 172 | if (fRegionList) delete [] fRegionList; |
89125dc7 | 173 | if (fMatList) delete fMatList; |
174 | if (fMatNames) {fMatNames->Delete(); delete fMatNames;} | |
8495a208 | 175 | if (gGeoManager) delete gGeoManager; |
176 | } | |
177 | ||
178 | // | |
179 | // private methods | |
180 | // | |
181 | //_____________________________________________________________________________ | |
182 | TFlukaMCGeometry::TFlukaMCGeometry(const TFlukaMCGeometry &) | |
89125dc7 | 183 | : TNamed() |
8495a208 | 184 | { |
185 | // | |
186 | // Copy constructor | |
187 | // | |
188 | } | |
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())) { | |
249 | imedium = vol->GetMedium()->GetId(); | |
250 | if (imedium == imed) { | |
251 | ireg = vol->GetNumber(); | |
252 | fRegionList[nreg++] = ireg; | |
253 | } | |
254 | } | |
255 | return fRegionList; | |
256 | } | |
257 | ||
258 | //_____________________________________________________________________________ | |
259 | Int_t *TFlukaMCGeometry::GetMaterialList(Int_t imat, Int_t &nreg) | |
260 | { | |
261 | // Get an ordered list of regions matching a given medium number | |
262 | nreg = 0; | |
263 | if (!fRegionList) fRegionList = new Int_t[NofVolumes()+1]; | |
264 | TIter next(gGeoManager->GetListOfUVolumes()); | |
265 | TGeoVolume *vol; | |
266 | Int_t imaterial, ireg; | |
267 | while ((vol = (TGeoVolume*)next())) { | |
268 | imaterial = vol->GetMedium()->GetMaterial()->GetIndex(); | |
269 | if (imaterial == imat) { | |
270 | ireg = vol->GetNumber(); | |
271 | fRegionList[nreg++] = ireg; | |
272 | } | |
273 | } | |
274 | return fRegionList; | |
275 | } | |
8495a208 | 276 | |
8495a208 | 277 | //_____________________________________________________________________________ |
278 | Int_t TFlukaMCGeometry::NofVolumes() const | |
279 | { | |
280 | // | |
281 | // Return total number of volumes in the geometry | |
282 | // | |
283 | ||
284 | return gGeoManager->GetListOfUVolumes()->GetEntriesFast()-1; | |
285 | } | |
8495a208 | 286 | |
287 | //_____________________________________________________________________________ | |
89125dc7 | 288 | TGeoMaterial * TFlukaMCGeometry::GetMakeWrongMaterial(Double_t z) |
289 | { | |
290 | // Try to replace a wrongly-defined material | |
291 | static Double_t kz[23] = {7.3, 17.8184, 7.2167, 10.856, 8.875, 8.9, 7.177, | |
292 | 25.72, 6.2363, 7.1315, 47.7056, 10.6467, 7.8598, 2.10853, 10.6001, 9.1193, | |
293 | 15.3383, 4.55, 9.6502, 6.4561, 21.7963, 29.8246, 15.4021}; | |
294 | ||
295 | Int_t ind; | |
296 | Double_t dz; | |
297 | for (ind=0; ind<23; ind++) { | |
298 | dz = TMath::Abs(z-kz[ind]); | |
299 | if (dz<1E-4) break; | |
300 | } | |
301 | if (ind>22) { | |
302 | printf("Cannot patch material with Z=%g\n", z); | |
303 | return 0; | |
304 | } | |
305 | TGeoMixture *mix = 0; | |
306 | TGeoElement *element; | |
81f1d030 | 307 | TGeoElementTable *table = gGeoManager->GetElementTable(); |
89125dc7 | 308 | switch (ind) { |
309 | case 0: // AIR | |
81f1d030 | 310 | mix = new TGeoMixture("TPC_AIR", 4, 0.001205); |
89125dc7 | 311 | element = table->GetElement(6); // C |
312 | mix->DefineElement(0, element, 0.000124); | |
313 | element = table->GetElement(7); // N | |
314 | mix->DefineElement(1, element, 0.755267); | |
315 | element = table->GetElement(8); // O | |
316 | mix->DefineElement(2, element, 0.231781); | |
317 | element = table->GetElement(18); // AR | |
318 | mix->DefineElement(3, element, 0.012827); | |
319 | break; | |
320 | case 1: //SDD SI CHIP | |
81f1d030 | 321 | mix = new TGeoMixture("ITS_SDD_SI", 6, 2.4485); |
89125dc7 | 322 | element = table->GetElement(1); |
323 | mix->DefineElement(0, element, 0.004367771); | |
324 | element = table->GetElement(6); | |
325 | mix->DefineElement(1, element, 0.039730642); | |
326 | element = table->GetElement(7); | |
327 | mix->DefineElement(2, element, 0.001396798); | |
328 | element = table->GetElement(8); | |
329 | mix->DefineElement(3, element, 0.01169634); | |
330 | element = table->GetElement(14); | |
331 | mix->DefineElement(4, element, 0.844665); | |
332 | element = table->GetElement(47); | |
333 | mix->DefineElement(5, element, 0.09814344903); | |
334 | break; | |
335 | case 2: // WATER | |
81f1d030 | 336 | mix = new TGeoMixture("ITS_WATER", 2, 1.0); |
89125dc7 | 337 | element = table->GetElement(1); |
338 | mix->DefineElement(0, element, 0.111898344); | |
339 | element = table->GetElement(8); | |
340 | mix->DefineElement(1, element, 0.888101656); | |
341 | break; | |
342 | case 3: // CERAMICS | |
81f1d030 | 343 | mix = new TGeoMixture("ITS_CERAMICS", 5, 3.6); |
89125dc7 | 344 | element = table->GetElement(8); |
345 | mix->DefineElement(0, element, 0.59956); | |
346 | element = table->GetElement(13); | |
347 | mix->DefineElement(1, element, 0.3776); | |
348 | element = table->GetElement(14); | |
349 | mix->DefineElement(2, element, 0.00933); | |
350 | element = table->GetElement(24); | |
351 | mix->DefineElement(3, element, 0.002); | |
352 | element = table->GetElement(25); | |
353 | mix->DefineElement(4, element, 0.0115); | |
354 | break; | |
355 | case 4: // EPOXY | |
81f1d030 | 356 | mix = new TGeoMixture("MUON_G10FR4", 4, 1.8); |
89125dc7 | 357 | element = table->GetElement(1); |
358 | mix->DefineElement(0, element, 0.19); | |
359 | element = table->GetElement(6); | |
360 | mix->DefineElement(1, element, 0.18); | |
361 | element = table->GetElement(8); | |
362 | mix->DefineElement(2, element, 0.35); | |
363 | element = table->GetElement(14); | |
364 | mix->DefineElement(3, element, 0.28); | |
365 | break; | |
366 | case 5: // EPOXY | |
367 | mix = new TGeoMixture("G10FR4", 4, 1.8); | |
368 | element = table->GetElement(1); | |
369 | mix->DefineElement(0, element, 0.19); | |
370 | element = table->GetElement(6); | |
371 | mix->DefineElement(1, element, 0.18); | |
372 | element = table->GetElement(8); | |
373 | mix->DefineElement(2, element, 0.35); | |
374 | element = table->GetElement(14); | |
375 | mix->DefineElement(3, element, 0.28); | |
376 | break; | |
377 | case 6: // KAPTON | |
81f1d030 | 378 | mix = new TGeoMixture("ITS_KAPTON", 4, 1.3); |
89125dc7 | 379 | element = table->GetElement(1); |
380 | mix->DefineElement(0, element, 0.026363415); | |
381 | element = table->GetElement(6); | |
382 | mix->DefineElement(1, element, 0.6911272); | |
383 | element = table->GetElement(7); | |
384 | mix->DefineElement(2, element, 0.073271325); | |
385 | element = table->GetElement(8); | |
386 | mix->DefineElement(3, element, 0.209238060); | |
387 | break; | |
388 | case 7: // INOX | |
81f1d030 | 389 | mix = new TGeoMixture("ITS_INOX", 9, 7.9); |
89125dc7 | 390 | element = table->GetElement(6); |
391 | mix->DefineElement(0, element, 0.0003); | |
392 | element = table->GetElement(14); | |
393 | mix->DefineElement(1, element, 0.01); | |
394 | element = table->GetElement(15); | |
395 | mix->DefineElement(2, element, 0.00045); | |
396 | element = table->GetElement(16); | |
397 | mix->DefineElement(3, element, 0.0003); | |
398 | element = table->GetElement(24); | |
399 | mix->DefineElement(4, element, 0.17); | |
400 | element = table->GetElement(25); | |
401 | mix->DefineElement(5, element, 0.02); | |
402 | element = table->GetElement(26); | |
403 | mix->DefineElement(6, element, 0.654); | |
404 | element = table->GetElement(28); | |
405 | mix->DefineElement(7, element, 0.12); | |
406 | element = table->GetElement(42); | |
407 | mix->DefineElement(8, element, 0.025); | |
408 | break; | |
409 | case 8: // ROHACELL | |
410 | mix = new TGeoMixture("ROHACELL", 4, 0.05); | |
411 | element = table->GetElement(1); | |
412 | mix->DefineElement(0, element, 0.07836617); | |
413 | element = table->GetElement(6); | |
414 | mix->DefineElement(1, element, 0.64648941); | |
415 | element = table->GetElement(7); | |
416 | mix->DefineElement(2, element, 0.08376983); | |
417 | element = table->GetElement(8); | |
418 | mix->DefineElement(3, element, 0.19137459); | |
419 | break; | |
420 | case 9: // SDD-C-AL | |
81f1d030 | 421 | mix = new TGeoMixture("ITS_SDD-C-AL", 5, 1.9837); |
89125dc7 | 422 | element = table->GetElement(1); |
423 | mix->DefineElement(0, element, 0.022632); | |
424 | element = table->GetElement(6); | |
425 | mix->DefineElement(1, element, 0.8176579); | |
426 | element = table->GetElement(7); | |
427 | mix->DefineElement(2, element, 0.0093488); | |
428 | element = table->GetElement(8); | |
429 | mix->DefineElement(3, element, 0.0503618); | |
430 | element = table->GetElement(13); | |
431 | mix->DefineElement(4, element, 0.1); | |
432 | break; | |
433 | case 10: // X7R-CAP | |
81f1d030 | 434 | mix = new TGeoMixture("ITS_X7R-CAP", 7, 6.72); |
89125dc7 | 435 | element = table->GetElement(8); |
436 | mix->DefineElement(0, element, 0.085975822); | |
437 | element = table->GetElement(22); | |
438 | mix->DefineElement(1, element, 0.084755042); | |
439 | element = table->GetElement(28); | |
440 | mix->DefineElement(2, element, 0.038244751); | |
441 | element = table->GetElement(29); | |
442 | mix->DefineElement(3, element, 0.009471271); | |
443 | element = table->GetElement(50); | |
444 | mix->DefineElement(4, element, 0.321736471); | |
445 | element = table->GetElement(56); | |
446 | mix->DefineElement(5, element, 0.251639432); | |
447 | element = table->GetElement(82); | |
448 | mix->DefineElement(6, element, 0.2081768); | |
449 | break; | |
450 | case 11: // SDD ruby sph. Al2O3 | |
81f1d030 | 451 | mix = new TGeoMixture("ITS_AL2O3", 2, 3.97); |
89125dc7 | 452 | element = table->GetElement(8); |
453 | mix->DefineElement(0, element, 0.5293); | |
454 | element = table->GetElement(13); | |
455 | mix->DefineElement(1, element, 0.4707); | |
456 | break; | |
457 | case 12: // SDD HV microcable | |
81f1d030 | 458 | mix = new TGeoMixture("ITS_HV-CABLE", 5, 1.6087); |
89125dc7 | 459 | element = table->GetElement(1); |
460 | mix->DefineElement(0, element, 0.01983871336); | |
461 | element = table->GetElement(6); | |
462 | mix->DefineElement(1, element, 0.520088819984); | |
463 | element = table->GetElement(7); | |
464 | mix->DefineElement(2, element, 0.0551367996); | |
465 | element = table->GetElement(8); | |
466 | mix->DefineElement(3, element, 0.157399667056); | |
467 | element = table->GetElement(13); | |
468 | mix->DefineElement(4, element, 0.247536); | |
469 | break; | |
470 | case 13: //SDD LV+signal cable | |
81f1d030 | 471 | mix = new TGeoMixture("ITS_LV-CABLE", 5, 2.1035); |
89125dc7 | 472 | element = table->GetElement(1); |
473 | mix->DefineElement(0, element, 0.0082859922); | |
474 | element = table->GetElement(6); | |
475 | mix->DefineElement(1, element, 0.21722436468); | |
476 | element = table->GetElement(7); | |
477 | mix->DefineElement(2, element, 0.023028867); | |
478 | element = table->GetElement(8); | |
479 | mix->DefineElement(3, element, 0.06574077612); | |
480 | element = table->GetElement(13); | |
481 | mix->DefineElement(4, element, 0.68572); | |
482 | break; | |
483 | case 14: //SDD hybrid microcab | |
81f1d030 | 484 | mix = new TGeoMixture("ITS_HYB-CAB", 5, 2.0502); |
89125dc7 | 485 | element = table->GetElement(1); |
486 | mix->DefineElement(0, element, 0.00926228815); | |
487 | element = table->GetElement(6); | |
488 | mix->DefineElement(1, element, 0.24281879711); | |
489 | element = table->GetElement(7); | |
490 | mix->DefineElement(2, element, 0.02574224025); | |
491 | element = table->GetElement(8); | |
492 | mix->DefineElement(3, element, 0.07348667449); | |
493 | element = table->GetElement(13); | |
494 | mix->DefineElement(4, element, 0.64869); | |
495 | break; | |
496 | case 15: //SDD anode microcab | |
81f1d030 | 497 | mix = new TGeoMixture("ITS_ANOD-CAB", 5, 1.7854); |
89125dc7 | 498 | element = table->GetElement(1); |
499 | mix->DefineElement(0, element, 0.0128595919215); | |
500 | element = table->GetElement(6); | |
501 | mix->DefineElement(1, element, 0.392653705471); | |
502 | element = table->GetElement(7); | |
503 | mix->DefineElement(2, element, 0.041626868025); | |
504 | element = table->GetElement(8); | |
505 | mix->DefineElement(3, element, 0.118832707289); | |
506 | element = table->GetElement(13); | |
507 | mix->DefineElement(4, element, 0.431909); | |
508 | break; | |
509 | case 16: // inox/alum | |
81f1d030 | 510 | mix = new TGeoMixture("ITS_INOX-AL", 5, 3.0705); |
89125dc7 | 511 | element = table->GetElement(13); |
512 | mix->DefineElement(0, element, 0.816164); | |
513 | element = table->GetElement(14); | |
514 | mix->DefineElement(1, element, 0.000919182); | |
515 | element = table->GetElement(24); | |
516 | mix->DefineElement(2, element, 0.0330906); | |
517 | element = table->GetElement(26); | |
518 | mix->DefineElement(3, element, 0.131443); | |
519 | element = table->GetElement(28); | |
520 | mix->DefineElement(4, element, 0.0183836); | |
521 | case 17: // MYLAR | |
81f1d030 | 522 | mix = new TGeoMixture("TPC_MYLAR", 3, 1.39); |
89125dc7 | 523 | element = table->GetElement(1); |
524 | mix->DefineElement(0, element, 0.0416667); | |
525 | element = table->GetElement(6); | |
526 | mix->DefineElement(1, element, 0.625); | |
527 | element = table->GetElement(8); | |
528 | mix->DefineElement(2, element, 0.333333); | |
529 | break; | |
530 | case 18: // SPDBUS(AL+KPT+EPOX) - unknown composition | |
81f1d030 | 531 | mix = new TGeoMixture("ITS_SPDBUS", 1, 1.906); |
89125dc7 | 532 | element = table->GetElement(9); |
533 | mix->DefineElement(0, element, 1.); | |
95161cd7 | 534 | z = element->Z(); |
89125dc7 | 535 | break; |
536 | case 19: // SDD/SSD rings - unknown composition | |
81f1d030 | 537 | mix = new TGeoMixture("ITS_SDDRINGS", 1, 1.8097); |
89125dc7 | 538 | element = table->GetElement(6); |
539 | mix->DefineElement(0, element, 1.); | |
95161cd7 | 540 | z = element->Z(); |
89125dc7 | 541 | break; |
542 | case 20: // SPD end ladder - unknown composition | |
81f1d030 | 543 | mix = new TGeoMixture("ITS_SPDEL", 1, 3.6374); |
89125dc7 | 544 | element = table->GetElement(22); |
545 | mix->DefineElement(0, element, 1.); | |
95161cd7 | 546 | z = element->Z(); |
89125dc7 | 547 | break; |
548 | case 21: // SDD end ladder - unknown composition | |
81f1d030 | 549 | mix = new TGeoMixture("ITS_SDDEL", 1, 0.3824); |
89125dc7 | 550 | element = table->GetElement(30); |
551 | mix->DefineElement(0, element, 1.); | |
95161cd7 | 552 | z = element->Z(); |
89125dc7 | 553 | break; |
554 | case 22: // SSD end ladder - unknown composition | |
81f1d030 | 555 | mix = new TGeoMixture("ITS_SSDEL", 1, 0.68); |
89125dc7 | 556 | element = table->GetElement(16); |
557 | mix->DefineElement(0, element, 1.); | |
95161cd7 | 558 | z = element->Z(); |
89125dc7 | 559 | break; |
560 | } | |
561 | mix->SetZ(z); | |
562 | printf("Patched with mixture %s\n", mix->GetName()); | |
563 | return mix; | |
8495a208 | 564 | } |
565 | ||
8495a208 | 566 | //_____________________________________________________________________________ |
efde9b4d | 567 | void TFlukaMCGeometry::CreateFlukaMatFile(const char *fname) |
8495a208 | 568 | { |
569 | // ==== from FLUGG ==== | |
570 | // NAMES OF ELEMENTS AND COMPOUNDS: the names must be written in upper case, | |
571 | // according to the fluka standard. In addition,. they must be equal to the | |
572 | // names of the fluka materials - see fluka manual - in order that the | |
573 | // program load the right cross sections, and equal to the names included in | |
574 | // the .pemf. Otherwise the user must define the LOW-MAT CARDS, and make his | |
575 | // own .pemf, in order to get the right cross sections loaded in memory. | |
89125dc7 | 576 | |
577 | ||
8495a208 | 578 | TString sname; |
579 | gGeoManager->Export("flgeom.root"); | |
580 | if (fname) sname = fname; | |
581 | else sname = "flukaMat.inp"; | |
582 | ofstream out; | |
583 | out.open(sname.Data(), ios::out); | |
584 | if (!out.good()) { | |
585 | Fatal("CreateFlukaMatFile", "could not open file %s for writing", sname.Data()); | |
586 | return; | |
587 | } | |
588 | PrintHeader(out, "MATERIALS AND COMPOUNDS"); | |
589 | PrintHeader(out, "MATERIALS"); | |
89125dc7 | 590 | Int_t i,j,idmat; |
591 | Int_t counttothree, nelem; | |
592 | Double_t a,z,rho, w; | |
81f1d030 | 593 | TGeoElementTable *table = gGeoManager->GetElementTable(); |
89125dc7 | 594 | TGeoElement *element; |
595 | element = table->GetElement(13); | |
596 | element->SetTitle("ALUMINUM"); // this is how FLUKA likes it ... | |
597 | element = table->GetElement(15); | |
598 | element->SetTitle("PHOSPHO"); // same story ... | |
599 | // element = table->GetElement(10); | |
600 | // element->SetTitle("ARGON"); // NEON not in neutron xsec table | |
601 | Int_t nelements = table->GetNelements(); | |
8495a208 | 602 | TList *matlist = gGeoManager->GetListOfMaterials(); |
95161cd7 | 603 | // TList *medlist = gGeoManager->GetListOfMedia(); |
604 | // Int_t nmed = medlist->GetSize(); | |
8495a208 | 605 | TIter next(matlist); |
606 | Int_t nmater = matlist->GetSize(); | |
607 | Int_t nfmater = 0; | |
89125dc7 | 608 | TGeoMaterial *mat; |
8495a208 | 609 | TGeoMixture *mix = 0; |
610 | TString matname; | |
89125dc7 | 611 | TObjString *objstr; |
612 | // Create all needed elements | |
613 | for (Int_t i=1; i<nelements; i++) { | |
614 | element = table->GetElement(i); | |
615 | // skip elements which are not defined | |
616 | if (!element->IsUsed() && !element->IsDefined()) continue; | |
617 | matname = element->GetTitle(); | |
618 | ToFlukaString(matname); | |
619 | rho = 0.999; | |
3860051a | 620 | |
89125dc7 | 621 | mat = new TGeoMaterial(matname, element->A(), element->Z(), rho); |
622 | mat->SetIndex(nfmater+3); | |
623 | mat->SetUsed(kTRUE); | |
624 | fMatList->Add(mat); | |
625 | objstr = new TObjString(matname.Data()); | |
626 | fMatNames->Add(objstr); | |
627 | nfmater++; | |
628 | } | |
6d184c54 | 629 | |
630 | fIndmat = nfmater; | |
95161cd7 | 631 | // TGeoMedium *med; |
89125dc7 | 632 | // Adjust material names and add them to FLUKA list |
8495a208 | 633 | for (i=0; i<nmater; i++) { |
634 | mat = (TGeoMaterial*)matlist->At(i); | |
81f1d030 | 635 | if (!mat->IsUsed()) continue; |
89125dc7 | 636 | z = mat->GetZ(); |
637 | a = mat->GetA(); | |
638 | rho = mat->GetDensity(); | |
639 | if (mat->GetZ()<0.001) { | |
efde9b4d | 640 | mat->SetIndex(2); // vacuum, built-in inside FLUKA |
641 | continue; | |
89125dc7 | 642 | } |
643 | matname = mat->GetName(); | |
644 | FlukaMatName(matname); | |
81f1d030 | 645 | |
89125dc7 | 646 | mat->SetIndex(nfmater+3); |
647 | objstr = new TObjString(matname.Data()); | |
648 | fMatList->Add(mat); | |
649 | fMatNames->Add(objstr); | |
650 | nfmater++; | |
651 | } | |
652 | ||
653 | // Dump all elements with MATERIAL cards | |
8495a208 | 654 | for (i=0; i<nfmater; i++) { |
89125dc7 | 655 | mat = (TGeoMaterial*)fMatList->At(i); |
656 | // mat->SetUsed(kFALSE); | |
657 | mix = 0; | |
8495a208 | 658 | out << setw(10) << "MATERIAL "; |
659 | out.setf(static_cast<std::ios::fmtflags>(0),std::ios::floatfield); | |
89125dc7 | 660 | objstr = (TObjString*)fMatNames->At(i); |
b0853588 | 661 | matname = objstr->GetString(); |
89125dc7 | 662 | z = mat->GetZ(); |
663 | a = mat->GetA(); | |
664 | rho = mat->GetDensity(); | |
8495a208 | 665 | if (mat->IsMixture()) { |
666 | out << setw(10) << " "; | |
667 | out << setw(10) << " "; | |
668 | mix = (TGeoMixture*)mat; | |
669 | } else { | |
89125dc7 | 670 | out << setw(10) << setiosflags(ios::fixed) << setprecision(1) << z; |
671 | out << setw(10) << setprecision(3) << a; | |
8495a208 | 672 | } |
673 | out.setf(static_cast<std::ios::fmtflags>(0),std::ios::floatfield); | |
89125dc7 | 674 | out << setw(10) << setiosflags(ios::scientific) << setprecision(3) << rho; |
8495a208 | 675 | out.setf(static_cast<std::ios::fmtflags>(0),std::ios::floatfield); |
89125dc7 | 676 | out << setw(10) << setiosflags(ios::fixed) << setprecision(1) << Double_t(mat->GetIndex()); |
8495a208 | 677 | out << setw(10) << " "; |
efde9b4d | 678 | out << setw(10) << " "; |
679 | out << setw(8) << matname.Data() << endl; | |
89125dc7 | 680 | if (!mix) { |
681 | // add LOW-MAT card for NEON to associate with ARGON neutron xsec | |
95161cd7 | 682 | if (z==10) { |
89125dc7 | 683 | out << setw(10) << "LOW-MAT "; |
684 | out.setf(static_cast<std::ios::fmtflags>(0),std::ios::floatfield); | |
685 | out << setw(10) << setiosflags(ios::fixed) << setprecision(1) << Double_t(mat->GetIndex()); | |
686 | out << setw(10) << setiosflags(ios::fixed) << setprecision(1) << 18.; | |
687 | out << setw(10) << setiosflags(ios::fixed) << setprecision(1) << -2.; | |
688 | out << setw(10) << setiosflags(ios::fixed) << setprecision(1) << 293.; | |
689 | out << setw(10) << " "; | |
690 | out << setw(10) << " "; | |
691 | // out << setw(8) << matname.Data() << endl; | |
692 | out << setw(8) << " " << endl; | |
95161cd7 | 693 | } |
694 | else { | |
695 | element = table->GetElement((int)z); | |
696 | TString elename = element->GetTitle(); | |
697 | ToFlukaString(elename); | |
698 | if( matname.CompareTo( elename ) != 0 ) { | |
699 | out << setw(10) << "LOW-MAT "; | |
700 | out.setf(static_cast<std::ios::fmtflags>(0),std::ios::floatfield); | |
701 | out << setw(10) << setiosflags(ios::fixed) << setprecision(1) << Double_t(mat->GetIndex()); | |
702 | out << setw(10) << setiosflags(ios::fixed) << setprecision(1) << z; | |
703 | out << setw(10) << setiosflags(ios::fixed) << setprecision(1) << " "; | |
704 | out << setw(10) << setiosflags(ios::fixed) << setprecision(1) << " "; | |
705 | out << setw(10) << " "; | |
706 | out << setw(10) << " "; | |
707 | // missing material at Low Energy Cross Section Table | |
708 | if( (int)z==10 || (int)z==21 || (int)z==34 || (int)z==37 || (int)z==39 || (int)z==44 || | |
709 | (int)z==45 || (int)z==46 || (int)z==52 || (int)z==57 || (int)z==59 || (int)z==60 || | |
710 | (int)z==61 || (int)z==65 || (int)z==66 || (int)z==67 || (int)z==68 || (int)z==69 || | |
711 | (int)z==70 || (int)z==71 || (int)z==72 || (int)z==76 || (int)z==77 || (int)z==78 || | |
712 | (int)z==81 || (int)z==84 || (int)z==85 || (int)z==86 || (int)z==87 || (int)z==88 || | |
713 | (int)z==89 || (int)z==91 ) | |
714 | out << setw(8) << "UNKNOWN " << endl; | |
715 | else | |
716 | out << setw(8) << elename.Data() << endl; | |
717 | // out << setw(8) << " " << endl; | |
718 | } | |
719 | } | |
89125dc7 | 720 | continue; |
d11bf3ca | 721 | } |
8495a208 | 722 | counttothree = 0; |
723 | out << setw(10) << "COMPOUND "; | |
724 | nelem = mix->GetNelements(); | |
89125dc7 | 725 | objstr = (TObjString*)fMatNames->At(i); |
8495a208 | 726 | matname = objstr->GetString(); |
8495a208 | 727 | for (j=0; j<nelem; j++) { |
89125dc7 | 728 | w = (mix->GetWmixt())[j]; |
729 | if (w<0.00001) w=0.00001; | |
730 | z = (mix->GetZmixt())[j]; | |
731 | a = (mix->GetAmixt())[j]; | |
732 | idmat = GetElementIndex(Int_t(z)); | |
733 | if (!idmat) Error("CreateFlukaMatFile", "element with Z=%f not found", z); | |
8495a208 | 734 | out.setf(static_cast<std::ios::fmtflags>(0),std::ios::floatfield); |
89125dc7 | 735 | out << setw(10) << setiosflags(ios::fixed) << setprecision(6) << -w; |
8495a208 | 736 | out.setf(static_cast<std::ios::fmtflags>(0),std::ios::floatfield); |
737 | out << setw(10) << setiosflags(ios::fixed) << setprecision(1) << Double_t(idmat); | |
738 | counttothree++; | |
739 | if (counttothree == 3) { | |
740 | out << matname.Data(); | |
741 | out << endl; | |
742 | if ( (j+1) != nelem) out << setw(10) << "COMPOUND "; | |
743 | counttothree = 0; | |
89125dc7 | 744 | } |
745 | } | |
8495a208 | 746 | if (nelem%3) { |
747 | for (j=0; j<(3-(nelem%3)); j++) | |
748 | out << setw(10) << " " << setw(10) << " "; | |
749 | out << matname.Data(); | |
750 | out << endl; | |
89125dc7 | 751 | } |
752 | } | |
8495a208 | 753 | Int_t nvols = gGeoManager->GetListOfUVolumes()->GetEntriesFast()-1; |
754 | TGeoVolume *vol; | |
8495a208 | 755 | // Now print the material assignments |
d11bf3ca | 756 | Double_t flagfield = 0.; |
757 | printf("#############################################################\n"); | |
89125dc7 | 758 | if (gFluka->IsFieldEnabled()) { |
d11bf3ca | 759 | flagfield = 1.; |
760 | printf("Magnetic field enabled\n"); | |
761 | } else printf("Magnetic field disabled\n"); | |
762 | printf("#############################################################\n"); | |
763 | ||
8495a208 | 764 | PrintHeader(out, "TGEO MATERIAL ASSIGNMENTS"); |
765 | for (i=1; i<=nvols; i++) { | |
3860051a | 766 | |
8495a208 | 767 | vol = gGeoManager->GetVolume(i); |
768 | mat = vol->GetMedium()->GetMaterial(); | |
769 | idmat = mat->GetIndex(); | |
89125dc7 | 770 | for (Int_t j=0; j<nfmater; j++) { |
771 | mat = (TGeoMaterial*)fMatList->At(j); | |
772 | if (mat->GetIndex() == idmat) mat->SetUsed(kTRUE); | |
773 | } | |
3860051a | 774 | |
775 | Float_t hasfield = (vol->GetMedium()->GetParam(1) > 0) ? flagfield : 0.; | |
89aa6d28 | 776 | out << "* Assigning material: " << vol->GetMedium()->GetMaterial()->GetName() << " to Volume: " << vol->GetName(); |
777 | out << endl; | |
778 | ||
8495a208 | 779 | out << setw(10) << "ASSIGNMAT "; |
780 | out.setf(static_cast<std::ios::fmtflags>(0),std::ios::floatfield); | |
781 | out << setw(10) << setiosflags(ios::fixed) << Double_t(idmat); | |
782 | out << setw(10) << setiosflags(ios::fixed) << Double_t(i); | |
783 | out << setw(10) << "0.0"; | |
b0853588 | 784 | out << setw(10) << "0.0"; |
3860051a | 785 | out << setw(10) << setiosflags(ios::fixed) << hasfield; |
b0853588 | 786 | out << setw(10) << "0.0"; |
8495a208 | 787 | out << endl; |
788 | } | |
8495a208 | 789 | out.close(); |
efde9b4d | 790 | fLastMaterial = nfmater+2; |
6d184c54 | 791 | } |
792 | ||
793 | void TFlukaMCGeometry::CreatePemfFile() | |
794 | { | |
795 | // | |
796 | // Steering routine to write and process peg files producing the pemf input | |
797 | // | |
798 | char number[20]; | |
799 | Int_t countMatOK = 0; | |
800 | Int_t countElemError = 0; | |
801 | Int_t countNoStern = 0; | |
802 | Int_t countMixError = 0; | |
803 | Int_t countGas = 0; | |
804 | Int_t countPemfError = 0; | |
805 | Int_t i; | |
806 | TGeoMaterial* mat = 0x0; | |
807 | TString sname; | |
808 | ||
809 | for (i = fIndmat; i < fLastMaterial - 2; i++) { | |
810 | printf("Write Peg Files %d\n", i); | |
811 | ||
812 | mat = (TGeoMaterial*)fMatList->At(i); | |
813 | if (!mat->IsUsed()) continue; | |
814 | sname = "mat"; | |
815 | sprintf(number, "%d", i); | |
816 | sname.Append(number); | |
817 | cout << endl; | |
818 | cout << endl; | |
819 | cout << "******************************************************************************" << endl; | |
820 | cout << "******************************************************************************" << endl; | |
821 | cout << endl; | |
822 | WritePegFile(i, &countNoStern, &countElemError, &countMixError, &countGas); | |
823 | sname.Prepend("$FLUPRO/pemf/rpemf peg/"); | |
824 | gSystem->Exec(sname.Data()); | |
825 | ||
826 | // check if the pemf file was created | |
827 | TString sname = Form("peg/mat%d.pemf", i); | |
828 | ifstream in( sname.Data() ); | |
829 | if ( in ) { | |
830 | countMatOK++; | |
831 | in.close(); | |
832 | } else { | |
833 | cout << "ERROR Fail to create the pemf file " << sname << endl; | |
834 | countPemfError++; | |
835 | } | |
836 | } | |
837 | cout << "Materials (pemf created) " << countMatOK << endl; | |
838 | cout << "Not Sternheimer par. found " << countNoStern << endl; | |
839 | cout << "Elements with error definitions (Z not integer) " << countElemError << endl; | |
840 | cout << "Mixtures with error definitions (Z not integer) " << countMixError << endl; | |
841 | cout << "Posible Gas (rho < 0.01) " << countGas << endl; | |
842 | // cout << "Posible Gas (without pressure information) " << countGasError << endl; | |
95161cd7 | 843 | cout << "Pemf files Error " << countPemfError << endl; |
6d184c54 | 844 | cout << endl << endl; |
845 | ||
846 | sname = "cat peg/*.pemf > peg/FlukaVmc.pemf"; | |
847 | gSystem->Exec(sname.Data()); | |
848 | sname = "mv peg/FlukaVmc.pemf FlukaVmc.pemf"; | |
849 | gSystem->Exec(sname.Data()); | |
89125dc7 | 850 | } |
851 | ||
852 | //_____________________________________________________________________________ | |
95161cd7 | 853 | void TFlukaMCGeometry::WritePegFile(Int_t imat, Int_t *NoStern, Int_t *ElemError, |
854 | Int_t *MixError, Int_t *countGas) const | |
89125dc7 | 855 | { |
95161cd7 | 856 | // Write the .peg file for one material |
857 | ||
89125dc7 | 858 | TGeoMaterial *mat = (TGeoMaterial*)fMatList->At(imat); |
859 | TString name = ((TObjString*)fMatNames->At(imat))->GetString(); | |
860 | TString line; | |
861 | char number[20]; | |
862 | TGeoElement *elem = mat->GetElement(); | |
863 | name = name.Strip(); | |
864 | TString sname = "mat"; | |
865 | sprintf(number, "%d", imat); | |
866 | sname.Append(number); | |
867 | sname.Append(".peg"); | |
868 | sname.Prepend("peg/"); | |
869 | ofstream out; | |
870 | out.open(sname.Data(), ios::out); | |
871 | if (!out.good()) return; | |
872 | Double_t dens = mat->GetDensity(); | |
873 | TGeoMixture *mix = 0; | |
874 | Int_t nel = 1; | |
875 | Int_t i; | |
876 | if (mat->IsMixture()) { | |
877 | mix = (TGeoMixture*)mat; | |
878 | nel = mix->GetNelements(); | |
95161cd7 | 879 | } |
880 | ||
89125dc7 | 881 | if (nel==1) { |
95161cd7 | 882 | cout << "( Element ) " << name << " Z=" << mat->GetZ() << " Rho " << mat->GetDensity() << endl; |
883 | ||
884 | Double_t zel = mat->GetZ(); | |
885 | if( (zel-Int_t(zel))>0.001 || zel < 1 ) { | |
886 | cout << " ERROR: A Element with not integer Z=" << zel << endl; | |
887 | cout << endl; | |
888 | (*ElemError)++; | |
889 | return; | |
890 | } | |
891 | ||
89125dc7 | 892 | out << "ELEM" << endl; |
893 | out << " &INP IRAYL=1, RHO=" << dens << ", " << endl; | |
95161cd7 | 894 | |
895 | // check for the Sternheimer parameters | |
896 | Double_t *issb_parm = GetISSB( mat->GetDensity(), 1, &zel, 0 ); | |
897 | if( issb_parm[0] > 0 && issb_parm[1] > 0 ) { | |
898 | cout << "Sternheimer parameters found" << endl; | |
899 | out << ", ISSB=1, IEV=" << issb_parm[0] << ", CBAR=" << issb_parm[1] | |
900 | << ", X0=" << issb_parm[2] << "," << endl; | |
901 | out << "X1=" <<issb_parm[3] <<", AFACT="<<issb_parm[4] <<", SK=" | |
902 | << issb_parm[5] << ", DELTA0=" << issb_parm[6]; | |
903 | } | |
904 | else { | |
905 | cout << "WARNING: Strange element, Sternheimer parameters not found" << endl; | |
906 | (*NoStern)++; | |
907 | } | |
908 | ||
909 | if (dens<0.01) { | |
910 | (*countGas)++; | |
911 | out << " GASP=1." << endl; | |
912 | } | |
913 | ||
89125dc7 | 914 | out << " &END" << endl; |
915 | out << name.Data() << endl; | |
916 | out << elem->GetName() << endl; | |
95161cd7 | 917 | |
918 | } | |
919 | else { | |
920 | ||
921 | cout << "( Mixture ) " << name << " Rho " << dens << " nElem " << nel << endl; | |
922 | ||
923 | Double_t *zt = new Double_t[nel]; | |
924 | Double_t *wt = new Double_t[nel]; | |
925 | for (int j=0; j<nel; j++) { | |
926 | zt[j] = (mix->GetZmixt())[j]; | |
927 | wt[j] = (mix->GetWmixt())[j]; | |
928 | if( (zt[j]-Int_t(zt[j])) > 0.001 || zt[j] < 1 ) { | |
929 | cout << "ERROR Mixture " << name << " with an element with not integer Z=" << zt[j] << endl; | |
930 | cout << endl; | |
931 | (*MixError)++; | |
932 | // just continue since the mixtures are not patch, | |
933 | // but the final release should include the return | |
934 | // return; | |
935 | } | |
936 | } | |
937 | Double_t *issb_parm = GetISSB( mat->GetDensity(), nel, zt, wt ); | |
89125dc7 | 938 | out << "MIXT" << endl; |
95161cd7 | 939 | out << " &INP IRAYL=1, NE=" << nel << ", RHOZ=" << wt[0] << ","; |
940 | line = Form(" &INP IRAYL=1, NE=%d RHOZ=%g", nel, wt[0]); | |
941 | for(int j=1; j<nel; j++) { | |
942 | out << " " << wt[j] << ","; | |
943 | line += Form(" %g,", wt[j] ); | |
944 | if( line.Length() > 60 ) { out << endl; line = ""; } | |
89125dc7 | 945 | } |
95161cd7 | 946 | out << " RHO=" << mat->GetDensity() << ", "; |
947 | line += Form(" RHO=%g, ", mat->GetDensity()); | |
948 | if( line.Length() > 60 ) { out << endl; line = ""; } | |
949 | ||
950 | if( issb_parm[0] > 0 && issb_parm[1] > 0 ) { | |
951 | cout << "Sternheimer parameters found" << endl; | |
952 | out << " ISSB=1, IEV=" << issb_parm[0] << ","; | |
953 | line += Form(" ISSB=1, IEV=%g,", issb_parm[0]); | |
954 | if( line.Length() > 60 ) { out << endl; line = ""; } | |
955 | out << " CBAR=" << issb_parm[1] << ","; | |
956 | line += Form(" CBAR=%g,",issb_parm[1]); | |
957 | if( line.Length() > 60 ) { out << endl; line = ""; } | |
958 | out << " X0=" << issb_parm[2] << ","; | |
959 | line += Form(" X0=%g,", issb_parm[2]); | |
960 | if( line.Length() > 60 ) { out << endl; line = ""; } | |
961 | out << " X1=" << issb_parm[3] << ","; | |
962 | line += Form(" X1=%g,", issb_parm[3]); | |
963 | if( line.Length() > 60 ) { out << endl; line = ""; } | |
964 | out << " AFACT="<< issb_parm[4] << ","; | |
965 | line += Form(" AFACT=%g,", issb_parm[4]); | |
966 | if( line.Length() > 60 ) { out << endl; line = ""; } | |
967 | out << " SK=" << issb_parm[5] << ","; | |
968 | line += Form(" SK=%g,", issb_parm[5]); | |
969 | if( line.Length() > 60 ) { out << endl; line = ""; } | |
970 | } | |
971 | else { | |
972 | cout << "Sternheimer parameters not found" << endl; | |
973 | (*NoStern)++; | |
974 | } | |
975 | ||
976 | if (dens<0.01){ | |
977 | (*countGas)++; | |
978 | out << " GASP=1." << endl; | |
979 | } | |
980 | ||
981 | out << " &END" << endl; | |
89125dc7 | 982 | out << name.Data() << endl; |
983 | for (i=0; i<nel; i++) { | |
984 | elem = mix->GetElement(i); | |
985 | line = elem->GetName(); | |
986 | if (line.Length()==1) line.Append(" "); | |
987 | out << line.Data() << " "; | |
988 | } | |
989 | out << endl; | |
95161cd7 | 990 | |
991 | delete [] zt; | |
992 | delete [] wt; | |
89125dc7 | 993 | } |
95161cd7 | 994 | |
6d184c54 | 995 | Double_t ue = 3000000.; // [MeV] |
996 | Double_t up = 3000000.; // [MeV] | |
997 | Double_t ae = -1.; | |
998 | Double_t ap = -1.; | |
999 | ||
1000 | ||
1001 | TObjArray* cutList = ((TFluka*) gMC)->GetListOfUserConfigs(); | |
1002 | TIter next(cutList); | |
1003 | TFlukaConfigOption* proc; | |
1004 | ||
1005 | while((proc = (TFlukaConfigOption*)next())) | |
1006 | { | |
1007 | if (proc->Medium() == mat->GetIndex()) { | |
1008 | ap = proc->Cut(kCUTGAM); | |
1009 | ae = proc->Cut(kCUTELE); | |
1010 | if (ap == -1.) ap = TFlukaConfigOption::DefaultCut(kCUTGAM); | |
1011 | if (ae == -1.) ae = TFlukaConfigOption::DefaultCut(kCUTELE); | |
1012 | break; | |
1013 | } | |
1014 | } | |
1015 | ||
1016 | if (ap == -1.) ap = TFlukaConfigOption::DefaultCut(kCUTGAM); | |
1017 | if (ae == -1.) ae = TFlukaConfigOption::DefaultCut(kCUTELE); | |
1018 | ||
1019 | ap *= 1000.; // [MeV] | |
1020 | ae = (ae + 0.00051099906) * 1000.; // [MeV] | |
1021 | ||
89125dc7 | 1022 | out << "ENER" << endl; |
6d184c54 | 1023 | out << " $INP AE=" << ae << ", UE=" << ue <<", AP=" << ap << ", UP=" << up << " $END" << endl; |
89125dc7 | 1024 | out << "PWLF" << endl; |
1025 | out << " $INP NALE=300, NALG=400, NALR=100 $END" << endl; | |
1026 | out << "DECK" << endl; | |
1027 | out << " $INP $END" << endl; | |
1028 | out << "TEST" << endl; | |
1029 | out << " $INP $END" << endl; | |
1030 | out.close(); | |
8495a208 | 1031 | } |
1032 | ||
95161cd7 | 1033 | Double_t * TFlukaMCGeometry::GetISSB(Double_t rho, Int_t nElem, Double_t *zelem, Double_t *welem ) const |
1034 | { | |
1035 | // Read the density effect parameters | |
1036 | // from R.M. Sternheimer et al. Atomic Data | |
1037 | // and Nuclear Data Tables, Vol. 30 No. 2 | |
1038 | // | |
1039 | // return the parameters if the element/mixture match with one of the list | |
1040 | // otherwise returns the parameters set to 0 | |
1041 | ||
1042 | struct sternheimerData { | |
1043 | TString longname; // element/mixture name | |
1044 | Int_t nelems; // number of constituents N | |
1045 | Int_t Z[20]; //[nelems] Z | |
1046 | Double_t wt[20]; //[nelems] weight fraction | |
1047 | Double_t density; // g/cm3 | |
1048 | Double_t iev; // Average Ion potential (eV) | |
6d184c54 | 1049 | // **** Sternheimer parameters **** |
95161cd7 | 1050 | Double_t cbar; // CBAR |
1051 | Double_t x0; // X0 | |
1052 | Double_t x1; // X1 | |
1053 | Double_t afact; // AFACT | |
1054 | Double_t sk; // SK | |
1055 | Double_t delta0; // DELTA0 | |
1056 | }; | |
1057 | ||
1058 | TString shortname; | |
1059 | TString formula; | |
1060 | Int_t num; | |
1061 | char state; | |
1062 | ||
1063 | static Double_t parameters[7]; | |
1064 | memset( parameters, 0, sizeof(Double_t) ); | |
1065 | ||
1066 | static sternheimerData sternDataArray[300]; | |
1067 | static Bool_t isFileRead = kFALSE; | |
1068 | ||
1069 | // Read the data file if is needed | |
1070 | if( isFileRead == kFALSE ) { | |
1071 | TString sSternheimerInp = getenv("ALICE_ROOT"); | |
1072 | sSternheimerInp +="/TFluka/input/Sternheimer.data"; | |
1073 | ||
1074 | ifstream in(sSternheimerInp); | |
1075 | char line[100]; | |
1076 | in.getline(line, 100); | |
1077 | in.getline(line, 100); | |
1078 | in.getline(line, 100); | |
1079 | in.getline(line, 100); | |
1080 | in.getline(line, 100); | |
1081 | in.getline(line, 100); | |
1082 | ||
1083 | ||
1084 | Int_t is = 0; | |
1085 | while( !in.eof() ) { | |
1086 | in >> shortname >> num >> sternDataArray[is].nelems | |
1087 | >> sternDataArray[is].longname >> formula >> state; | |
1088 | if( in.eof() ) break; | |
1089 | for(int i=0; i<sternDataArray[is].nelems; i++) { | |
1090 | in >> sternDataArray[is].Z[i] >> sternDataArray[is].wt[i]; | |
1091 | } | |
1092 | in >> sternDataArray[is].density; | |
1093 | in >> sternDataArray[is].iev; | |
1094 | in >> sternDataArray[is].cbar; | |
1095 | in >> sternDataArray[is].x0; | |
1096 | in >> sternDataArray[is].x1; | |
1097 | in >> sternDataArray[is].afact; | |
1098 | in >> sternDataArray[is].sk; | |
1099 | if( sternDataArray[is].nelems == 1 ) in >> sternDataArray[is].delta0; | |
1100 | is++; | |
1101 | } | |
1102 | isFileRead = kTRUE; | |
1103 | in.close(); | |
1104 | } | |
1105 | ||
1106 | Int_t is = 0; | |
1107 | while( is < 280 ) { | |
1108 | ||
1109 | // check for elements | |
1110 | if( sternDataArray[is].nelems == 1 && nElem == 1 | |
1111 | && sternDataArray[is].Z[0] == Int_t(*zelem) | |
1112 | && TMath::Abs( (sternDataArray[is].density - rho)/sternDataArray[is].density ) < 0.1 ) { | |
1113 | cout << sternDataArray[is].longname << " #elems:" << sternDataArray[is].nelems << " Rho:" | |
1114 | << sternDataArray[is].density << endl; | |
1115 | cout << sternDataArray[is].iev << " " | |
1116 | << sternDataArray[is].cbar << " " | |
1117 | << sternDataArray[is].x0 << " " | |
1118 | << sternDataArray[is].x1 << " " | |
1119 | << sternDataArray[is].afact << " " | |
1120 | << sternDataArray[is].sk << " " | |
1121 | << sternDataArray[is].delta0 << endl; | |
1122 | ||
1123 | parameters[0] = sternDataArray[is].iev; | |
1124 | parameters[1] = sternDataArray[is].cbar; | |
1125 | parameters[2] = sternDataArray[is].x0; | |
1126 | parameters[3] = sternDataArray[is].x1; | |
1127 | parameters[4] = sternDataArray[is].afact; | |
1128 | parameters[5] = sternDataArray[is].sk; | |
1129 | parameters[6] = sternDataArray[is].delta0; | |
1130 | return parameters; | |
1131 | } | |
1132 | ||
1133 | // check for mixture | |
1134 | int nmatch = 0; | |
1135 | if( sternDataArray[is].nelems > 1 && sternDataArray[is].nelems == nElem ) { | |
1136 | for(int j=0; j<sternDataArray[is].nelems; j++) { | |
1137 | if( sternDataArray[is].Z[j] == Int_t(zelem[j]) && | |
1138 | TMath::Abs( (sternDataArray[is].wt[j] - welem[j])/sternDataArray[is].wt[j] ) < 0.1 ) | |
1139 | nmatch++; | |
1140 | } | |
1141 | } | |
1142 | ||
1143 | if( sternDataArray[is].nelems > 1 && | |
1144 | TMath::Abs( (sternDataArray[is].density - rho)/sternDataArray[is].density ) < 0.1 | |
1145 | && nmatch == sternDataArray[is].nelems ) { | |
1146 | cout << sternDataArray[is].longname << " #elem:" << sternDataArray[is].nelems << " Rho:" | |
1147 | << sternDataArray[is].density << endl; | |
1148 | cout << sternDataArray[is].iev << " " | |
1149 | << sternDataArray[is].cbar << " " | |
1150 | << sternDataArray[is].x0 << " " | |
1151 | << sternDataArray[is].x1 << " " | |
1152 | << sternDataArray[is].afact << " " | |
1153 | << sternDataArray[is].sk << " " | |
1154 | << sternDataArray[is].delta0 << endl; | |
1155 | ||
1156 | parameters[0] = sternDataArray[is].iev; | |
1157 | parameters[1] = sternDataArray[is].cbar; | |
1158 | parameters[2] = sternDataArray[is].x0; | |
1159 | parameters[3] = sternDataArray[is].x1; | |
1160 | parameters[4] = sternDataArray[is].afact; | |
1161 | parameters[5] = sternDataArray[is].sk; | |
1162 | parameters[6] = 0; | |
1163 | return parameters; | |
1164 | } | |
1165 | is++; | |
1166 | } | |
1167 | return parameters; | |
1168 | } | |
1169 | ||
8495a208 | 1170 | //_____________________________________________________________________________ |
1171 | void TFlukaMCGeometry::PrintHeader(ofstream &out, const char *text) const | |
1172 | { | |
1173 | // Print a FLUKA header. | |
1174 | out << "*\n" << "*\n" << "*\n"; | |
1175 | out << "********************* " << text << " *********************\n" | |
1176 | << "*\n"; | |
1177 | out << "*...+....1....+....2....+....3....+....4....+....5....+....6....+....7..." | |
1178 | << endl; | |
1179 | out << "*" << endl; | |
1180 | } | |
1181 | ||
1182 | //_____________________________________________________________________________ | |
1183 | Int_t TFlukaMCGeometry::RegionId() const | |
1184 | { | |
1185 | // Returns current region id <-> TGeo node id | |
1186 | if (gGeoManager->IsOutside()) return 0; | |
1187 | return gGeoManager->GetCurrentNode()->GetUniqueID(); | |
1188 | } | |
89125dc7 | 1189 | |
1190 | //_____________________________________________________________________________ | |
1191 | Int_t TFlukaMCGeometry::GetElementIndex(Int_t z) const | |
1192 | { | |
6f5667d1 | 1193 | // Get index of a material having a given Z element. |
89125dc7 | 1194 | TIter next(fMatList); |
1195 | TGeoMaterial *mat; | |
1196 | Int_t index = 0; | |
1197 | while ((mat=(TGeoMaterial*)next())) { | |
1198 | if (mat->IsMixture()) continue; | |
1199 | if (mat->GetElement()->Z() == z) return mat->GetIndex(); | |
1200 | } | |
1201 | return index; | |
1202 | } | |
1203 | ||
05265ca9 | 1204 | //_____________________________________________________________________________ |
1205 | void TFlukaMCGeometry::SetMreg(Int_t mreg) | |
1206 | { | |
1207 | // Update if needed next history; | |
89125dc7 | 1208 | if (gFluka->GetDummyBoundary()==2) { |
d11bf3ca | 1209 | gGeoManager->CdNode(fNextLattice-1); |
1210 | return; | |
1211 | } | |
89125dc7 | 1212 | Int_t curreg = (gGeoManager->IsOutside())?(gMCGeom->NofVolumes()+1):gGeoManager->GetCurrentVolume()->GetNumber(); |
05265ca9 | 1213 | if (mreg==curreg) return; |
1214 | if (mreg==fNextRegion) { | |
1215 | if (fNextLattice!=999999999) gGeoManager->CdNode(fNextLattice-1); | |
1216 | return; | |
d11bf3ca | 1217 | } else { |
1218 | if (mreg == fCurrentRegion) { | |
1219 | if (fCurrentLattice!=999999999) gGeoManager->CdNode(fCurrentLattice-1); | |
1220 | return; | |
1221 | } | |
1222 | } | |
2bc4c610 | 1223 | if (fDebug) printf("ERROR: mreg=%i neither current nor next region\n", mreg); |
05265ca9 | 1224 | } |
1225 | ||
d11bf3ca | 1226 | //_____________________________________________________________________________ |
1227 | void TFlukaMCGeometry::SetCurrentRegion(Int_t mreg, Int_t latt) | |
1228 | { | |
1229 | // Set index/history for next entered region | |
1230 | fCurrentRegion = mreg; | |
1231 | fCurrentLattice = latt; | |
1232 | } | |
1233 | ||
05265ca9 | 1234 | //_____________________________________________________________________________ |
1235 | void TFlukaMCGeometry::SetNextRegion(Int_t mreg, Int_t latt) | |
1236 | { | |
1237 | // Set index/history for next entered region | |
1238 | fNextRegion = mreg; | |
1239 | fNextLattice = latt; | |
1240 | } | |
8495a208 | 1241 | |
1242 | //_____________________________________________________________________________ | |
1243 | void TFlukaMCGeometry::ToFlukaString(TString &str) const | |
1244 | { | |
1245 | // ToFlukaString converts an string to something usefull in FLUKA: | |
1246 | // * Capital letters | |
1247 | // * Only 8 letters | |
1248 | // * Replace ' ' by '_' | |
1249 | if (str.Length()<8) { | |
1250 | str += " "; | |
1251 | } | |
1252 | str.Remove(8); | |
1253 | Int_t ilast; | |
1254 | for (ilast=7; ilast>0; ilast--) if (str(ilast)!=' ') break; | |
1255 | str.ToUpper(); | |
1256 | for (Int_t pos=0; pos<ilast; pos++) | |
1257 | if (str(pos)==' ') str.Replace(pos,1,"_",1); | |
1258 | return; | |
1259 | } | |
89125dc7 | 1260 | |
1261 | //_____________________________________________________________________________ | |
1262 | void TFlukaMCGeometry::FlukaMatName(TString &str) const | |
1263 | { | |
81f1d030 | 1264 | // Strip the detector name |
1265 | ||
1266 | TObjArray * tokens = str.Tokenize("_"); | |
1267 | Int_t ntok = tokens->GetEntries(); | |
1268 | if (ntok > 0) { | |
1269 | TString head = ((TObjString*) tokens->At(0))->GetString(); | |
1270 | Int_t nhead = head.Length(); | |
1271 | str = str.Remove(0, nhead + 1); | |
1272 | } | |
1273 | tokens->Clear(); | |
1274 | delete tokens; | |
1275 | ||
6f5667d1 | 1276 | // Convert a name to upper case 8 chars. |
89125dc7 | 1277 | ToFlukaString(str); |
1278 | Int_t ilast; | |
1279 | for (ilast=7; ilast>0; ilast--) if (str(ilast)!=' ') break; | |
1280 | if (ilast>5) ilast = 5; | |
1281 | char number[3]; | |
1282 | TIter next(fMatNames); | |
1283 | TObjString *objstr; | |
1284 | TString matname; | |
1285 | Int_t index = 0; | |
1286 | while ((objstr=(TObjString*)next())) { | |
1287 | matname = objstr->GetString(); | |
1288 | if (matname == str) { | |
1289 | index++; | |
1290 | if (index<10) { | |
1291 | number[0] = '0'; | |
1292 | sprintf(&number[1], "%d", index); | |
1293 | } else if (index<100) { | |
1294 | sprintf(number, "%d", index); | |
1295 | } else { | |
1296 | Error("FlukaMatName", "Too many materials %s", str.Data()); | |
1297 | return; | |
1298 | } | |
1299 | str.Replace(ilast+1, 2, number); | |
1300 | str.Remove(8); | |
1301 | } | |
1302 | } | |
1303 | } | |
1304 | ||
8495a208 | 1305 | //______________________________________________________________________________ |
1306 | void TFlukaMCGeometry::Vname(const char *name, char *vname) const | |
1307 | { | |
1308 | // | |
1309 | // convert name to upper case. Make vname at least 4 chars | |
1310 | // | |
1311 | Int_t l = strlen(name); | |
1312 | Int_t i; | |
1313 | l = l < 4 ? l : 4; | |
1314 | for (i=0;i<l;i++) vname[i] = toupper(name[i]); | |
1315 | for (i=l;i<4;i++) vname[i] = ' '; | |
1316 | vname[4] = 0; | |
1317 | } | |
1318 | ||
1319 | ||
81f1d030 | 1320 | |
8495a208 | 1321 | // FLUKA GEOMETRY WRAPPERS - to replace FLUGG wrappers |
1322 | ||
8495a208 | 1323 | //_____________________________________________________________________________ |
1324 | Int_t idnrwr(const Int_t & /*nreg*/, const Int_t & /*mlat*/) | |
1325 | { | |
1326 | // from FLUGG: | |
1327 | // Wrapper for setting DNEAR option on fluka side. Must return 0 | |
1328 | // if user doesn't want Fluka to use DNEAR to compute the | |
1329 | // step (the same effect is obtained with the GLOBAL (WHAT(3)=-1) | |
1330 | // card in fluka input), returns 1 if user wants Fluka always to | |
1331 | // use DNEAR (in this case, be sure that GEANT4 DNEAR is unique, | |
1332 | // coming from all directions!!!) | |
89125dc7 | 1333 | if (gMCGeom->IsDebugging()) printf("========== Dummy IDNRWR\n"); |
8495a208 | 1334 | return 0; |
1335 | } | |
1336 | ||
8495a208 | 1337 | //_____________________________________________________________________________ |
1338 | void g1wr(Double_t &pSx, Double_t &pSy, Double_t &pSz, | |
d11bf3ca | 1339 | Double_t *pV, Int_t &oldReg , const Int_t &oldLttc, Double_t &propStep, |
cefd9d9b | 1340 | Int_t &/*nascFlag*/, Double_t &retStep, Int_t &newReg, |
2573ac89 | 1341 | Double_t &saf, Int_t &newLttc, Int_t <tcFlag, |
8495a208 | 1342 | Double_t *sLt, Int_t *jrLt) |
d11bf3ca | 1343 | |
8495a208 | 1344 | { |
d11bf3ca | 1345 | // Initialize FLUKa point and direction; |
6f5667d1 | 1346 | gNstep++; |
784c3bdd | 1347 | NORLAT.xn[0] = pSx; |
1348 | NORLAT.xn[1] = pSy; | |
1349 | NORLAT.xn[2] = pSz; | |
1350 | NORLAT.wn[0] = pV[0]; | |
1351 | NORLAT.wn[1] = pV[1]; | |
1352 | NORLAT.wn[2] = pV[2]; | |
8495a208 | 1353 | |
89125dc7 | 1354 | if (gMCGeom->IsDebugging()) { |
2bc4c610 | 1355 | printf("========== Inside G1WR\n"); |
1356 | printf(" point/dir:(%14.9f, %14.9f, %14.9f, %g, %g, %g)\n", pSx,pSy,pSz,pV[0],pV[1],pV[2]); | |
d11bf3ca | 1357 | printf(" oldReg=%i oldLttc=%i pstep=%f\n",oldReg, oldLttc, propStep); |
a940d455 | 1358 | } |
1359 | Int_t olttc = oldLttc; | |
81f1d030 | 1360 | if (oldLttc<=0) { |
a940d455 | 1361 | gGeoManager->FindNode(pSx,pSy,pSz); |
1362 | olttc = gGeoManager->GetCurrentNodeId()+1; | |
81f1d030 | 1363 | if (gMCGeom->IsDebugging()) { |
1364 | printf("WOOPS: old reg/latt = %i/%i\n",oldReg,oldLttc); | |
1365 | printf("point: (%16.12f, %16.12f, %16.12f) in lttc=%i\n", pSx,pSy,pSz,olttc); | |
1366 | } | |
2bc4c610 | 1367 | } |
81f1d030 | 1368 | Int_t ccreg,cclat; |
1369 | gMCGeom->GetCurrentRegion(ccreg,cclat); | |
a940d455 | 1370 | gMCGeom->SetCurrentRegion(oldReg, olttc); |
d11bf3ca | 1371 | // Initialize default return values |
1372 | lttcFlag = 0; | |
a940d455 | 1373 | jrLt[lttcFlag] = olttc; |
d11bf3ca | 1374 | sLt[lttcFlag] = propStep; |
1375 | jrLt[lttcFlag+1] = -1; | |
1376 | sLt[lttcFlag+1] = 0.; | |
1377 | newReg = oldReg; | |
a940d455 | 1378 | newLttc = olttc; |
d11bf3ca | 1379 | // check if dummy boundary flag is set |
1380 | Int_t curLttc, curReg; | |
89125dc7 | 1381 | if (gFluka->IsDummyBoundary()) { |
cefd9d9b | 1382 | // printf("Dummy boundary intercepted. Point is: %f, %f, %f\n", pSx, pSy, pSz); |
a940d455 | 1383 | Bool_t crossedDummy = (olttc == TFlukaMCGeometry::kLttcVirtual)?kTRUE:kFALSE; |
d11bf3ca | 1384 | if (crossedDummy) { |
1385 | // FLUKA crossed the dummy boundary - update new region/history | |
1386 | retStep = 0.; | |
1387 | saf = 0.; | |
89125dc7 | 1388 | gMCGeom->GetNextRegion(newReg, newLttc); |
1389 | gMCGeom->SetMreg(newReg); | |
1390 | if (gMCGeom->IsDebugging()) printf(" virtual newReg=%i newLttc=%i\n", newReg, newLttc); | |
d11bf3ca | 1391 | sLt[lttcFlag] = 0.; // null step in current region |
1392 | lttcFlag++; | |
1393 | jrLt[lttcFlag] = newLttc; | |
1394 | sLt[lttcFlag] = 0.; // null step in next region | |
1395 | jrLt[lttcFlag+1] = -1; | |
1396 | sLt[lttcFlag+1] = 0.; | |
89125dc7 | 1397 | gFluka->SetDummyBoundary(0); |
d11bf3ca | 1398 | return; |
1399 | } | |
1400 | } | |
1401 | ||
1402 | // Reset outside flag | |
cefd9d9b | 1403 | gGeoManager->SetOutside(kFALSE); |
d11bf3ca | 1404 | |
1405 | // Reset dummy boundary flag | |
89125dc7 | 1406 | gFluka->SetDummyBoundary(0); |
d11bf3ca | 1407 | |
1408 | curLttc = gGeoManager->GetCurrentNodeId()+1; | |
1409 | curReg = gGeoManager->GetCurrentVolume()->GetNumber(); | |
a940d455 | 1410 | if (olttc != curLttc) { |
d11bf3ca | 1411 | // FLUKA crossed the boundary : we trust that the given point is really there, |
1412 | // so we just update TGeo state | |
a940d455 | 1413 | gGeoManager->CdNode(olttc-1); |
2573ac89 | 1414 | curLttc = gGeoManager->GetCurrentNodeId()+1; |
d11bf3ca | 1415 | curReg = gGeoManager->GetCurrentVolume()->GetNumber(); |
89125dc7 | 1416 | if (gMCGeom->IsDebugging()) printf(" re-initialized point: curReg=%i curLttc=%i\n", curReg, curLttc); |
d11bf3ca | 1417 | } |
784c3bdd | 1418 | // Now the current TGeo state reflects the FLUKA state |
1419 | Double_t extra = 1.E-10; | |
1420 | ||
1421 | ||
89125dc7 | 1422 | if (gMCGeom->IsDebugging()) printf(" current path: %s\n", gGeoManager->GetPath()); |
cefd9d9b | 1423 | gGeoManager->SetCurrentPoint(pSx+extra*pV[0], pSy+extra*pV[1], pSz+extra*pV[2]); |
1424 | gGeoManager->SetCurrentDirection(pV); | |
81f1d030 | 1425 | // gGeoManager->FindNextBoundary(-propStep); |
1426 | gGeoManager->FindNextBoundary(-100000.); | |
d11bf3ca | 1427 | Double_t snext = gGeoManager->GetStep(); |
784c3bdd | 1428 | |
1429 | if (snext<=0.0) { | |
81f1d030 | 1430 | snext = 0.0; |
1431 | /* | |
1432 | if (!crossed) { | |
1433 | // artefact due to MS | |
1434 | saf = 0.0; | |
1435 | newReg = -3; | |
1436 | sLt[lttcFlag] = 0.0; | |
1437 | if (gMCGeom->IsDebugging()) { | |
1438 | printf("BACK SCATTERING (reg/lttc = %i/%i)\n", oldReg,oldLttc); | |
1439 | printf("point: (%16.11f, %16.11f, %16.11f)\n", pSx,pSy,pSz); | |
1440 | printf("dir : (%16.11f, %16.11f, %16.11f)\n", pV[0],pV[1],pV[2]); | |
1441 | } | |
1442 | return; | |
1443 | } | |
1444 | */ | |
1445 | // else ... | |
1446 | // FLUKA native detects some extra errors even when a boundary was crossed, returning a -3 | |
1447 | // If the particle turns back to the original region before crossing at the first step, | |
1448 | // we just return the distance to the boundary, not issuing an error (not due to MS) | |
784c3bdd | 1449 | } |
cefd9d9b | 1450 | snext += extra; |
d11bf3ca | 1451 | saf = gGeoManager->GetSafeDistance(); |
cefd9d9b | 1452 | saf -= extra; |
1453 | if (saf<0) saf=0.0; | |
784c3bdd | 1454 | else saf -= saf*3.0e-09; |
81f1d030 | 1455 | // saf *= 0.3; |
1456 | PAREM.dist = snext; | |
784c3bdd | 1457 | NORLAT.distn = snext; |
1458 | NORLAT.xn[0] += snext*pV[0]; | |
1459 | NORLAT.xn[1] += snext*pV[1]; | |
1460 | NORLAT.xn[2] += snext*pV[2]; | |
cefd9d9b | 1461 | if (snext>propStep) { |
d11bf3ca | 1462 | // Next boundary further than proposed step, which is approved |
1463 | retStep = propStep; | |
1464 | sLt[lttcFlag] = propStep; | |
1465 | return; | |
1466 | } | |
1467 | // The next boundary is closer. We try to cross it. | |
81f1d030 | 1468 | if (saf>propStep) saf = propStep; // Safety should be less than the proposed step if a boundary will be crossed |
1469 | // saf = 0.0; // !!! SAFETY SHOULD BE 0 IF THE BOUNDARY WILL BE CROSSED ??? | |
cefd9d9b | 1470 | gGeoManager->SetCurrentPoint(pSx,pSy,pSz); |
2573ac89 | 1471 | Double_t *point = gGeoManager->GetCurrentPoint(); |
1472 | Double_t *dir = gGeoManager->GetCurrentDirection(); | |
d11bf3ca | 1473 | Double_t pt[3]; |
1474 | memcpy(pt, point, 3*sizeof(Double_t)); | |
1475 | ||
2573ac89 | 1476 | Int_t i; |
81f1d030 | 1477 | extra = 1.E-13; |
1478 | for (i=0;i<3;i++) point[i] += (snext+extra)*dir[i]; | |
cefd9d9b | 1479 | // locate next region |
d11bf3ca | 1480 | gGeoManager->FindNode(); |
1481 | newLttc = (gGeoManager->IsOutside())?(TFlukaMCGeometry::kLttcOutside):gGeoManager->GetCurrentNodeId()+1; | |
81f1d030 | 1482 | while (newLttc==olttc) { |
1483 | extra *= 10.; | |
1484 | if (extra>1.E-5) break; | |
1485 | for (i=0;i<3;i++) point[i] += extra*dir[i]; | |
1486 | gGeoManager->FindNode(); | |
1487 | newLttc = (gGeoManager->IsOutside())?(TFlukaMCGeometry::kLttcOutside):gGeoManager->GetCurrentNodeId()+1; | |
1488 | } | |
d11bf3ca | 1489 | gGeoManager->SetCurrentPoint(pt); |
89125dc7 | 1490 | newReg = (gGeoManager->IsOutside())?(gMCGeom->NofVolumes()+1):gGeoManager->GetCurrentVolume()->GetNumber(); |
1491 | if (gMCGeom->IsDebugging()) printf(" newReg=%i newLttc=%i\n", newReg, newLttc); | |
d11bf3ca | 1492 | |
81f1d030 | 1493 | |
d11bf3ca | 1494 | // We really crossed the boundary, but is it the same region ? |
89125dc7 | 1495 | gMCGeom->SetNextRegion(newReg, newLttc); |
a940d455 | 1496 | if (newReg==oldReg && newLttc!=olttc) { |
d11bf3ca | 1497 | // Virtual boundary between replicants |
89125dc7 | 1498 | if (gMCGeom->IsDebugging()) printf(" DUMMY boundary\n"); |
d11bf3ca | 1499 | newReg = 1; // cheat FLUKA telling it it crossed the TOP region |
1500 | newLttc = TFlukaMCGeometry::kLttcVirtual; | |
1501 | // mark that next boundary is virtual | |
89125dc7 | 1502 | gFluka->SetDummyBoundary(1); |
d11bf3ca | 1503 | } |
1504 | retStep = snext; | |
1505 | sLt[lttcFlag] = snext; | |
1506 | lttcFlag++; | |
1507 | jrLt[lttcFlag] = newLttc; | |
1508 | sLt[lttcFlag] = snext; | |
1509 | jrLt[lttcFlag+1] = -1; | |
1510 | sLt[lttcFlag+1] = 0.; | |
1511 | ||
a940d455 | 1512 | if (newLttc!=olttc) { |
cefd9d9b | 1513 | if (gGeoManager->IsOutside()) gGeoManager->SetOutside(kFALSE); |
a940d455 | 1514 | gGeoManager->CdNode(olttc-1); |
8495a208 | 1515 | } |
89125dc7 | 1516 | if (gMCGeom->IsDebugging()) { |
d11bf3ca | 1517 | printf("=> snext=%g safe=%g\n", snext, saf); |
1518 | for (Int_t i=0; i<lttcFlag+1; i++) printf(" jrLt[%i]=%i sLt[%i]=%g\n", i,jrLt[i],i,sLt[i]); | |
1519 | } | |
89125dc7 | 1520 | if (gMCGeom->IsDebugging()) printf("<= G1WR (in: %s)\n", gGeoManager->GetPath()); |
8495a208 | 1521 | } |
1522 | ||
efde9b4d | 1523 | //_____________________________________________________________________________ |
1524 | void g1rtwr() | |
1525 | { | |
89125dc7 | 1526 | if (gMCGeom->IsDebugging()) printf("========== Dummy G1RTWR\n"); |
efde9b4d | 1527 | } |
1528 | ||
1529 | //_____________________________________________________________________________ | |
1530 | void conhwr(Int_t & /*intHist*/, Int_t * /*incrCount*/) | |
1531 | { | |
89125dc7 | 1532 | if (gMCGeom->IsDebugging()) printf("========== Dummy CONHWR\n"); |
efde9b4d | 1533 | } |
1534 | ||
1535 | //_____________________________________________________________________________ | |
2573ac89 | 1536 | void inihwr(Int_t &intHist) |
efde9b4d | 1537 | { |
89125dc7 | 1538 | if (gMCGeom->IsDebugging()) printf("========== Inside INIHWR -> reinitializing history: %i\n", intHist); |
2573ac89 | 1539 | if (gGeoManager->IsOutside()) gGeoManager->CdTop(); |
1540 | if (intHist<=0) { | |
1541 | // printf("=== wrong history number\n"); | |
1542 | return; | |
1543 | } | |
1544 | if (intHist==0) gGeoManager->CdTop(); | |
1545 | else gGeoManager->CdNode(intHist-1); | |
89125dc7 | 1546 | if (gMCGeom->IsDebugging()) { |
2bc4c610 | 1547 | printf(" --- current path: %s\n", gGeoManager->GetPath()); |
1548 | printf("<= INIHWR\n"); | |
1549 | } | |
efde9b4d | 1550 | } |
1551 | ||
1552 | //_____________________________________________________________________________ | |
1553 | void jomiwr(const Int_t & /*nge*/, const Int_t & /*lin*/, const Int_t & /*lou*/, | |
1554 | Int_t &flukaReg) | |
1555 | { | |
1556 | // Geometry initialization wrapper called by FLUKAM. Provides to FLUKA the | |
1557 | // number of regions (volumes in TGeo) | |
1558 | // build application geometry | |
89125dc7 | 1559 | if (gMCGeom->IsDebugging()) printf("========== Inside JOMIWR\n"); |
2573ac89 | 1560 | flukaReg = gGeoManager->GetListOfUVolumes()->GetEntriesFast(); |
89125dc7 | 1561 | if (gMCGeom->IsDebugging()) printf("<= JOMIWR: last region=%i\n", flukaReg); |
efde9b4d | 1562 | } |
1563 | ||
1564 | //_____________________________________________________________________________ | |
1565 | void lkdbwr(Double_t &pSx, Double_t &pSy, Double_t &pSz, | |
784c3bdd | 1566 | Double_t *pV, const Int_t &oldReg, const Int_t &oldLttc, |
efde9b4d | 1567 | Int_t &newReg, Int_t &flagErr, Int_t &newLttc) |
1568 | { | |
89125dc7 | 1569 | if (gMCGeom->IsDebugging()) { |
2bc4c610 | 1570 | printf("========== Inside LKDBWR (%f, %f, %f)\n",pSx, pSy, pSz); |
784c3bdd | 1571 | printf(" in: pV=(%f, %f, %f)\n", pV[0], pV[1], pV[2]); |
2bc4c610 | 1572 | printf(" in: oldReg=%i oldLttc=%i\n", oldReg, oldLttc); |
1573 | } | |
784c3bdd | 1574 | lkwr(pSx,pSy,pSz,pV,oldReg,oldLttc,newReg,flagErr,newLttc); |
efde9b4d | 1575 | } |
1576 | ||
1577 | //_____________________________________________________________________________ | |
1578 | void lkfxwr(Double_t &pSx, Double_t &pSy, Double_t &pSz, | |
784c3bdd | 1579 | Double_t *pV, const Int_t &oldReg, const Int_t &oldLttc, |
efde9b4d | 1580 | Int_t &newReg, Int_t &flagErr, Int_t &newLttc) |
1581 | { | |
89125dc7 | 1582 | if (gMCGeom->IsDebugging()) { |
2bc4c610 | 1583 | printf("========== Inside LKFXWR (%f, %f, %f)\n",pSx, pSy, pSz); |
784c3bdd | 1584 | printf(" in: pV=(%f, %f, %f)\n", pV[0], pV[1], pV[2]); |
2bc4c610 | 1585 | printf(" in: oldReg=%i oldLttc=%i\n", oldReg, oldLttc); |
1586 | } | |
784c3bdd | 1587 | lkwr(pSx,pSy,pSz,pV,oldReg,oldLttc,newReg,flagErr,newLttc); |
efde9b4d | 1588 | } |
1589 | ||
1590 | //_____________________________________________________________________________ | |
1591 | void lkmgwr(Double_t &pSx, Double_t &pSy, Double_t &pSz, | |
784c3bdd | 1592 | Double_t *pV, const Int_t &oldReg, const Int_t &oldLttc, |
efde9b4d | 1593 | Int_t &flagErr, Int_t &newReg, Int_t &newLttc) |
1594 | { | |
89125dc7 | 1595 | if (gMCGeom->IsDebugging()) { |
2bc4c610 | 1596 | printf("========== Inside LKMGWR (%f, %f, %f)\n",pSx, pSy, pSz); |
784c3bdd | 1597 | printf(" in: pV=(%f, %f, %f)\n", pV[0], pV[1], pV[2]); |
2bc4c610 | 1598 | printf(" in: oldReg=%i oldLttc=%i\n", oldReg, oldLttc); |
1599 | } | |
784c3bdd | 1600 | lkwr(pSx,pSy,pSz,pV,oldReg,oldLttc,newReg,flagErr,newLttc); |
efde9b4d | 1601 | } |
1602 | ||
1603 | //_____________________________________________________________________________ | |
1604 | void lkwr(Double_t &pSx, Double_t &pSy, Double_t &pSz, | |
784c3bdd | 1605 | Double_t *pV, const Int_t &oldReg, const Int_t &oldLttc, |
efde9b4d | 1606 | Int_t &newReg, Int_t &flagErr, Int_t &newLttc) |
1607 | { | |
89125dc7 | 1608 | if (gMCGeom->IsDebugging()) { |
2bc4c610 | 1609 | printf("========== Inside LKWR (%f, %f, %f)\n",pSx, pSy, pSz); |
784c3bdd | 1610 | printf(" in: pV=(%f, %f, %f)\n", pV[0], pV[1], pV[2]); |
2bc4c610 | 1611 | printf(" in: oldReg=%i oldLttc=%i\n", oldReg, oldLttc); |
1612 | } | |
81f1d030 | 1613 | /* |
784c3bdd | 1614 | NORLAT.xn[0] = pSx; |
1615 | NORLAT.xn[1] = pSy; | |
1616 | NORLAT.xn[2] = pSz; | |
1617 | NORLAT.wn[0] = pV[0]; | |
1618 | NORLAT.wn[1] = pV[1]; | |
1619 | NORLAT.wn[2] = pV[2]; | |
81f1d030 | 1620 | */ |
efde9b4d | 1621 | TGeoNode *node = gGeoManager->FindNode(pSx, pSy, pSz); |
1622 | if (gGeoManager->IsOutside()) { | |
89125dc7 | 1623 | newReg = gMCGeom->NofVolumes()+1; |
05265ca9 | 1624 | // newLttc = gGeoManager->GetCurrentNodeId(); |
1625 | newLttc = 999999999; | |
89125dc7 | 1626 | if (gMCGeom->IsDebugging()) { |
2bc4c610 | 1627 | printf("OUTSIDE\n"); |
1628 | printf(" out: newReg=%i newLttc=%i\n", newReg, newLttc); | |
1629 | printf("<= LKMGWR\n"); | |
1630 | } | |
05265ca9 | 1631 | flagErr = newReg; |
1632 | return; | |
efde9b4d | 1633 | } |
1634 | newReg = node->GetVolume()->GetNumber(); | |
d11bf3ca | 1635 | newLttc = gGeoManager->GetCurrentNodeId()+1; |
89125dc7 | 1636 | gMCGeom->SetNextRegion(newReg, newLttc); |
efde9b4d | 1637 | flagErr = newReg; |
89125dc7 | 1638 | if (gMCGeom->IsDebugging()) { |
2bc4c610 | 1639 | printf(" out: newReg=%i newLttc=%i in %s\n", newReg, newLttc, gGeoManager->GetPath()); |
1640 | printf("<= LKWR\n"); | |
1641 | } | |
efde9b4d | 1642 | } |
1643 | ||
1644 | //_____________________________________________________________________________ | |
2573ac89 | 1645 | void nrmlwr(Double_t &pSx, Double_t &pSy, Double_t &pSz, |
1646 | Double_t &pVx, Double_t &pVy, Double_t &pVz, | |
1647 | Double_t *norml, const Int_t &oldReg, | |
1648 | const Int_t &newReg, Int_t &flagErr) | |
efde9b4d | 1649 | { |
89125dc7 | 1650 | if (gMCGeom->IsDebugging()) { |
2bc4c610 | 1651 | printf("========== Inside NRMLWR (%g, %g, %g, %g, %g, %g)\n", pSx,pSy,pSz,pVx,pVy,pVz); |
1652 | printf(" oldReg=%i, newReg=%i\n", oldReg,newReg); | |
1653 | } | |
784c3bdd | 1654 | gGeoManager->SetCurrentPoint(NORLAT.xn[0], NORLAT.xn[1], NORLAT.xn[2]); |
1655 | gGeoManager->SetCurrentDirection(NORLAT.wn[0], NORLAT.wn[1], NORLAT.wn[2]); | |
2573ac89 | 1656 | Double_t *dnorm = gGeoManager->FindNormalFast(); |
1657 | flagErr = 0; | |
1658 | if (!dnorm) { | |
1659 | printf(" ERROR: Cannot compute fast normal\n"); | |
1660 | flagErr = 1; | |
1661 | norml[0] = -pVx; | |
1662 | norml[1] = -pVy; | |
1663 | norml[2] = -pVz; | |
1664 | } | |
a668e34d | 1665 | else { |
1666 | norml[0] = -dnorm[0]; | |
1667 | norml[1] = -dnorm[1]; | |
1668 | norml[2] = -dnorm[2]; | |
1669 | } | |
89125dc7 | 1670 | if (gMCGeom->IsDebugging()) { |
784c3bdd | 1671 | printf(" normal to boundary: (%g, %g, %g)\n", norml[0], norml[1], norml[2]); |
2bc4c610 | 1672 | printf("<= NRMLWR\n"); |
1673 | } | |
efde9b4d | 1674 | } |
1675 | ||
1676 | //_____________________________________________________________________________ | |
1677 | void rgrpwr(const Int_t & /*flukaReg*/, const Int_t & /*ptrLttc*/, Int_t & /*g4Reg*/, | |
1678 | Int_t * /*indMother*/, Int_t * /*repMother*/, Int_t & /*depthFluka*/) | |
1679 | { | |
89125dc7 | 1680 | if (gMCGeom->IsDebugging()) printf("=> Dummy RGRPWR\n"); |
efde9b4d | 1681 | } |
1682 | ||
1683 | //_____________________________________________________________________________ | |
1684 | Int_t isvhwr(const Int_t &check, const Int_t & intHist) | |
1685 | { | |
1686 | // from FLUGG: | |
1687 | // Wrapper for saving current navigation history (fCheck=default) | |
1688 | // and returning its pointer. If fCheck=-1 copy of history pointed | |
1689 | // by intHist is made in NavHistWithCount object, and its pointer | |
1690 | // is returned. fCheck=1 and fCheck=2 cases are only in debugging | |
1691 | // version: an array is created by means of FGeometryInit functions | |
1692 | // (but could be a static int * ptrArray = new int[10000] with | |
1693 | // file scope as well) that stores a flag for deleted/undeleted | |
1694 | // histories and at the end of event is checked to verify that | |
1695 | // all saved history objects have been deleted. | |
1696 | ||
1697 | // For TGeo, just return the current node ID. No copy need to be made. | |
1698 | ||
89125dc7 | 1699 | if (gMCGeom->IsDebugging()) printf("=> Inside ISVHWR\n"); |
efde9b4d | 1700 | if (check<0) return intHist; |
2573ac89 | 1701 | Int_t histInt = gGeoManager->GetCurrentNodeId()+1; |
89125dc7 | 1702 | if (gMCGeom->IsDebugging()) printf("<= ISVHWR: history is: %i in: %s\n", histInt, gGeoManager->GetPath()); |
efde9b4d | 1703 | return histInt; |
1704 | } | |
1705 | ||
1706 | ||
8495a208 | 1707 | |
1708 |