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