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