1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
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
24 // Author: Andrei Gheata 10/07/2003
26 #include "Riostream.h"
31 #include "TFlukaMCGeometry.h"
32 #include "TFlukaConfigOption.h"
33 #include "TGeoManager.h"
34 #include "TGeoVolume.h"
35 #include "TObjString.h"
38 #include "Fstepsz.h" //(STEPSZ) fluka common
41 # define idnrwr idnrwr_
43 # define g1rtwr g1rtwr_
44 # define conhwr conhwr_
45 # define inihwr inihwr_
46 # define jomiwr jomiwr_
47 # define lkdbwr lkdbwr_
48 # define lkfxwr lkfxwr_
49 # define lkmgwr lkmgwr_
51 # define magfld magfld_
52 # define nrmlwr nrmlwr_
53 # define rgrpwr rgrpwr_
54 # define isvhwr isvhwr_
58 # define idnrwr IDNRWR
60 # define g1rtwr G1RTWR
61 # define conhwr CONHWR
62 # define inihwr INIHWR
63 # define jomiwr JOMIWR
64 # define lkdbwr LKDBWR
65 # define lkfxwr LKFXWR
66 # define lkmgwr LKMGWR
68 # define magfld MAGFLD
69 # define nrmlwr NRMLWR
70 # define rgrpwr RGRPWR
71 # define isvhwr ISVHWR
75 //____________________________________________________________________________
79 // Prototypes for FLUKA navigation methods
81 Int_t type_of_call idnrwr(const Int_t & /*nreg*/, const Int_t & /*mlat*/);
82 void type_of_call g1wr(Double_t & /*pSx*/, Double_t & /*pSy*/, Double_t & /*pSz*/,
83 Double_t * /*pV*/, Int_t & /*oldReg*/ , const Int_t & /*oldLttc*/, Double_t & /*propStep*/,
84 Int_t & /*nascFlag*/, Double_t & /*retStep*/, Int_t & /*newReg*/,
85 Double_t & /*saf*/, Int_t & /*newLttc*/, Int_t & /*LttcFlag*/,
86 Double_t *s /*Lt*/, Int_t * /*jrLt*/);
88 void type_of_call g1rtwr();
89 void type_of_call conhwr(Int_t & /*intHist*/, Int_t & /*incrCount*/);
90 void type_of_call inihwr(Int_t & /*intHist*/);
91 void type_of_call jomiwr(const Int_t & /*nge*/, const Int_t & /*lin*/, const Int_t & /*lou*/,
92 Int_t & /*flukaReg*/);
93 void type_of_call lkdbwr(Double_t & /*pSx*/, Double_t & /*pSy*/, Double_t & /*pSz*/,
94 Double_t * /*pV*/, const Int_t & /*oldReg*/, const Int_t & /*oldLttc*/,
95 Int_t & /*flagErr*/, Int_t & /*newReg*/, Int_t & /*newLttc*/);
96 void type_of_call lkfxwr(Double_t & /*pSx*/, Double_t & /*pSy*/, Double_t & /*pSz*/,
97 Double_t * /*pV*/, const Int_t & /*oldReg*/, const Int_t & /*oldLttc*/,
98 Int_t & /*flagErr*/, Int_t & /*newReg*/, Int_t & /*newLttc*/);
99 void type_of_call lkmgwr(Double_t & /*pSx*/, Double_t & /*pSy*/, Double_t & /*pSz*/,
100 Double_t * /*pV*/, const Int_t & /*oldReg*/, const Int_t & /*oldLttc*/,
101 Int_t & /*flagErr*/, Int_t & /*newReg*/, Int_t & /*newLttc*/);
102 void type_of_call lkwr(Double_t & /*pSx*/, Double_t & /*pSy*/, Double_t & /*pSz*/,
103 Double_t * /*pV*/, const Int_t & /*oldReg*/, const Int_t & /*oldLttc*/,
104 Int_t & /*flagErr*/, Int_t & /*newReg*/, Int_t & /*newLttc*/);
105 // void type_of_call magfld(const Double_t & /*pX*/, const Double_t & /*pY*/, const Double_t & /*pZ*/,
106 // Double_t & /*cosBx*/, Double_t & /*cosBy*/, Double_t & /*cosBz*/,
107 // Double_t & /*Bmag*/, Int_t & /*reg*/, Int_t & /*idiscflag*/);
108 void type_of_call nrmlwr(Double_t & /*pSx*/, Double_t & /*pSy*/, Double_t & /*pSz*/,
109 Double_t & /*pVx*/, Double_t & /*pVy*/, Double_t & /*pVz*/,
110 Double_t * /*norml*/, const Int_t & /*oldReg*/,
111 const Int_t & /*newReg*/, Int_t & /*flagErr*/);
112 void type_of_call rgrpwr(const Int_t & /*flukaReg*/, const Int_t & /*ptrLttc*/, Int_t & /*g4Reg*/,
113 Int_t * /*indMother*/, Int_t * /*repMother*/, Int_t & /*depthFluka*/);
114 Int_t type_of_call isvhwr(const Int_t & /*fCheck*/, const Int_t & /*intHist*/);
117 // TFluka global pointer
119 TFlukaMCGeometry *gMCGeom = 0;
122 ClassImp(TFlukaMCGeometry)
124 TFlukaMCGeometry* TFlukaMCGeometry::fgInstance= NULL;
126 //_____________________________________________________________________________
127 TFlukaMCGeometry::TFlukaMCGeometry(const char *name, const char *title)
128 :TNamed(name, title),
138 fMatList(new TObjArray(256)),
139 fMatNames(new TObjArray(256))
142 // Standard constructor
144 gFluka = (TFluka*)gMC;
149 //_____________________________________________________________________________
150 TFlukaMCGeometry::TFlukaMCGeometry()
166 // Default constructor
168 gFluka = (TFluka*)gMC;
173 //_____________________________________________________________________________
174 TFlukaMCGeometry::~TFlukaMCGeometry()
180 if (fRegionList) delete [] fRegionList;
181 if (fMatList) delete fMatList;
182 if (fMatNames) {fMatNames->Delete(); delete fMatNames;}
183 if (gGeoManager) delete gGeoManager;
190 //_____________________________________________________________________________
191 Double_t* TFlukaMCGeometry::CreateDoubleArray(Float_t* array, Int_t size) const
193 // Converts Float_t* array to Double_t*,
194 // !! The new array has to be deleted by user.
197 Double_t* doubleArray;
199 doubleArray = new Double_t[size];
200 for (Int_t i=0; i<size; i++) doubleArray[i] = array[i];
204 doubleArray = new Double_t[1];
212 //_____________________________________________________________________________
213 Int_t TFlukaMCGeometry::GetMedium() const
215 // Get current medium number
217 TGeoNode *node = gGeoManager->GetCurrentNode();
218 if (!node) imed = gGeoManager->GetTopNode()->GetVolume()->GetMedium()->GetId();
219 else imed = node->GetVolume()->GetMedium()->GetId();
220 if (fDebug) printf("GetMedium=%i\n", imed);
224 //_____________________________________________________________________________
225 Int_t TFlukaMCGeometry::GetFlukaMaterial(Int_t imed) const
227 // Returns FLUKA material index for medium IMED
228 TGeoMedium *med = (TGeoMedium*)gGeoManager->GetListOfMedia()->At(imed-1);
230 Error("GetFlukaMaterial", "MEDIUM %i nor found", imed);
233 TGeoMaterial* mat = med->GetMaterial();
234 if (!mat->IsUsed()) return -1;
235 Int_t imatfl = med->GetMaterial()->GetIndex();
239 //_____________________________________________________________________________
240 Int_t *TFlukaMCGeometry::GetRegionList(Int_t imed, Int_t &nreg)
242 // Get an ordered list of regions matching a given medium number
244 if (!fRegionList) fRegionList = new Int_t[NofVolumes()+1];
245 TIter next(gGeoManager->GetListOfUVolumes());
248 while ((vol = (TGeoVolume*)next())) {
250 if ((med = vol->GetMedium()) == 0) continue;
251 imedium = med->GetId();
252 if (imedium == imed) {
253 ireg = vol->GetNumber();
254 fRegionList[nreg++] = ireg;
260 //_____________________________________________________________________________
261 Int_t *TFlukaMCGeometry::GetMaterialList(Int_t imat, Int_t &nreg)
263 // Get an ordered list of regions matching a given medium number
265 if (!fRegionList) fRegionList = new Int_t[NofVolumes()+1];
266 TIter next(gGeoManager->GetListOfUVolumes());
268 Int_t imaterial, ireg;
269 while ((vol = (TGeoVolume*)next())) {
271 if ((med = vol->GetMedium()) == 0) continue;
272 imaterial = med->GetMaterial()->GetIndex();
273 if (imaterial == imat) {
274 ireg = vol->GetNumber();
275 fRegionList[nreg++] = ireg;
281 //_____________________________________________________________________________
282 Int_t TFlukaMCGeometry::NofVolumes() const
285 // Return total number of volumes in the geometry
288 return gGeoManager->GetListOfUVolumes()->GetEntriesFast()-1;
291 //_____________________________________________________________________________
292 TGeoMaterial * TFlukaMCGeometry::GetMakeWrongMaterial(Double_t z)
294 // Try to replace a wrongly-defined material
295 static Double_t kz[23] = {7.3, 17.8184, 7.2167, 10.856, 8.875, 8.9, 7.177,
296 25.72, 6.2363, 7.1315, 47.7056, 10.6467, 7.8598, 2.10853, 10.6001, 9.1193,
297 15.3383, 4.55, 9.6502, 6.4561, 21.7963, 29.8246, 15.4021};
301 for (ind=0; ind<23; ind++) {
302 dz = TMath::Abs(z-kz[ind]);
306 printf("Cannot patch material with Z=%g\n", z);
309 TGeoMixture *mix = 0;
310 TGeoElement *element;
311 TGeoElementTable *table = gGeoManager->GetElementTable();
314 mix = new TGeoMixture("TPC_AIR", 4, 0.001205);
315 element = table->GetElement(6); // C
316 mix->DefineElement(0, element, 0.000124);
317 element = table->GetElement(7); // N
318 mix->DefineElement(1, element, 0.755267);
319 element = table->GetElement(8); // O
320 mix->DefineElement(2, element, 0.231781);
321 element = table->GetElement(18); // AR
322 mix->DefineElement(3, element, 0.012827);
324 case 1: //SDD SI CHIP
325 mix = new TGeoMixture("ITS_SDD_SI", 6, 2.4485);
326 element = table->GetElement(1);
327 mix->DefineElement(0, element, 0.004367771);
328 element = table->GetElement(6);
329 mix->DefineElement(1, element, 0.039730642);
330 element = table->GetElement(7);
331 mix->DefineElement(2, element, 0.001396798);
332 element = table->GetElement(8);
333 mix->DefineElement(3, element, 0.01169634);
334 element = table->GetElement(14);
335 mix->DefineElement(4, element, 0.844665);
336 element = table->GetElement(47);
337 mix->DefineElement(5, element, 0.09814344903);
340 mix = new TGeoMixture("ITS_WATER", 2, 1.0);
341 element = table->GetElement(1);
342 mix->DefineElement(0, element, 0.111898344);
343 element = table->GetElement(8);
344 mix->DefineElement(1, element, 0.888101656);
347 mix = new TGeoMixture("ITS_CERAMICS", 5, 3.6);
348 element = table->GetElement(8);
349 mix->DefineElement(0, element, 0.59956);
350 element = table->GetElement(13);
351 mix->DefineElement(1, element, 0.3776);
352 element = table->GetElement(14);
353 mix->DefineElement(2, element, 0.00933);
354 element = table->GetElement(24);
355 mix->DefineElement(3, element, 0.002);
356 element = table->GetElement(25);
357 mix->DefineElement(4, element, 0.0115);
360 mix = new TGeoMixture("MUON_G10FR4", 4, 1.8);
361 element = table->GetElement(1);
362 mix->DefineElement(0, element, 0.19);
363 element = table->GetElement(6);
364 mix->DefineElement(1, element, 0.18);
365 element = table->GetElement(8);
366 mix->DefineElement(2, element, 0.35);
367 element = table->GetElement(14);
368 mix->DefineElement(3, element, 0.28);
371 mix = new TGeoMixture("G10FR4", 4, 1.8);
372 element = table->GetElement(1);
373 mix->DefineElement(0, element, 0.19);
374 element = table->GetElement(6);
375 mix->DefineElement(1, element, 0.18);
376 element = table->GetElement(8);
377 mix->DefineElement(2, element, 0.35);
378 element = table->GetElement(14);
379 mix->DefineElement(3, element, 0.28);
382 mix = new TGeoMixture("ITS_KAPTON", 4, 1.3);
383 element = table->GetElement(1);
384 mix->DefineElement(0, element, 0.026363415);
385 element = table->GetElement(6);
386 mix->DefineElement(1, element, 0.6911272);
387 element = table->GetElement(7);
388 mix->DefineElement(2, element, 0.073271325);
389 element = table->GetElement(8);
390 mix->DefineElement(3, element, 0.209238060);
393 mix = new TGeoMixture("ITS_INOX", 9, 7.9);
394 element = table->GetElement(6);
395 mix->DefineElement(0, element, 0.0003);
396 element = table->GetElement(14);
397 mix->DefineElement(1, element, 0.01);
398 element = table->GetElement(15);
399 mix->DefineElement(2, element, 0.00045);
400 element = table->GetElement(16);
401 mix->DefineElement(3, element, 0.0003);
402 element = table->GetElement(24);
403 mix->DefineElement(4, element, 0.17);
404 element = table->GetElement(25);
405 mix->DefineElement(5, element, 0.02);
406 element = table->GetElement(26);
407 mix->DefineElement(6, element, 0.654);
408 element = table->GetElement(28);
409 mix->DefineElement(7, element, 0.12);
410 element = table->GetElement(42);
411 mix->DefineElement(8, element, 0.025);
414 mix = new TGeoMixture("ROHACELL", 4, 0.05);
415 element = table->GetElement(1);
416 mix->DefineElement(0, element, 0.07836617);
417 element = table->GetElement(6);
418 mix->DefineElement(1, element, 0.64648941);
419 element = table->GetElement(7);
420 mix->DefineElement(2, element, 0.08376983);
421 element = table->GetElement(8);
422 mix->DefineElement(3, element, 0.19137459);
425 mix = new TGeoMixture("ITS_SDD-C-AL", 5, 1.9837);
426 element = table->GetElement(1);
427 mix->DefineElement(0, element, 0.022632);
428 element = table->GetElement(6);
429 mix->DefineElement(1, element, 0.8176579);
430 element = table->GetElement(7);
431 mix->DefineElement(2, element, 0.0093488);
432 element = table->GetElement(8);
433 mix->DefineElement(3, element, 0.0503618);
434 element = table->GetElement(13);
435 mix->DefineElement(4, element, 0.1);
438 mix = new TGeoMixture("ITS_X7R-CAP", 7, 6.72);
439 element = table->GetElement(8);
440 mix->DefineElement(0, element, 0.085975822);
441 element = table->GetElement(22);
442 mix->DefineElement(1, element, 0.084755042);
443 element = table->GetElement(28);
444 mix->DefineElement(2, element, 0.038244751);
445 element = table->GetElement(29);
446 mix->DefineElement(3, element, 0.009471271);
447 element = table->GetElement(50);
448 mix->DefineElement(4, element, 0.321736471);
449 element = table->GetElement(56);
450 mix->DefineElement(5, element, 0.251639432);
451 element = table->GetElement(82);
452 mix->DefineElement(6, element, 0.2081768);
454 case 11: // SDD ruby sph. Al2O3
455 mix = new TGeoMixture("ITS_AL2O3", 2, 3.97);
456 element = table->GetElement(8);
457 mix->DefineElement(0, element, 0.5293);
458 element = table->GetElement(13);
459 mix->DefineElement(1, element, 0.4707);
461 case 12: // SDD HV microcable
462 mix = new TGeoMixture("ITS_HV-CABLE", 5, 1.6087);
463 element = table->GetElement(1);
464 mix->DefineElement(0, element, 0.01983871336);
465 element = table->GetElement(6);
466 mix->DefineElement(1, element, 0.520088819984);
467 element = table->GetElement(7);
468 mix->DefineElement(2, element, 0.0551367996);
469 element = table->GetElement(8);
470 mix->DefineElement(3, element, 0.157399667056);
471 element = table->GetElement(13);
472 mix->DefineElement(4, element, 0.247536);
474 case 13: //SDD LV+signal cable
475 mix = new TGeoMixture("ITS_LV-CABLE", 5, 2.1035);
476 element = table->GetElement(1);
477 mix->DefineElement(0, element, 0.0082859922);
478 element = table->GetElement(6);
479 mix->DefineElement(1, element, 0.21722436468);
480 element = table->GetElement(7);
481 mix->DefineElement(2, element, 0.023028867);
482 element = table->GetElement(8);
483 mix->DefineElement(3, element, 0.06574077612);
484 element = table->GetElement(13);
485 mix->DefineElement(4, element, 0.68572);
487 case 14: //SDD hybrid microcab
488 mix = new TGeoMixture("ITS_HYB-CAB", 5, 2.0502);
489 element = table->GetElement(1);
490 mix->DefineElement(0, element, 0.00926228815);
491 element = table->GetElement(6);
492 mix->DefineElement(1, element, 0.24281879711);
493 element = table->GetElement(7);
494 mix->DefineElement(2, element, 0.02574224025);
495 element = table->GetElement(8);
496 mix->DefineElement(3, element, 0.07348667449);
497 element = table->GetElement(13);
498 mix->DefineElement(4, element, 0.64869);
500 case 15: //SDD anode microcab
501 mix = new TGeoMixture("ITS_ANOD-CAB", 5, 1.7854);
502 element = table->GetElement(1);
503 mix->DefineElement(0, element, 0.0128595919215);
504 element = table->GetElement(6);
505 mix->DefineElement(1, element, 0.392653705471);
506 element = table->GetElement(7);
507 mix->DefineElement(2, element, 0.041626868025);
508 element = table->GetElement(8);
509 mix->DefineElement(3, element, 0.118832707289);
510 element = table->GetElement(13);
511 mix->DefineElement(4, element, 0.431909);
513 case 16: // inox/alum
514 mix = new TGeoMixture("ITS_INOX-AL", 5, 3.0705);
515 element = table->GetElement(13);
516 mix->DefineElement(0, element, 0.816164);
517 element = table->GetElement(14);
518 mix->DefineElement(1, element, 0.000919182);
519 element = table->GetElement(24);
520 mix->DefineElement(2, element, 0.0330906);
521 element = table->GetElement(26);
522 mix->DefineElement(3, element, 0.131443);
523 element = table->GetElement(28);
524 mix->DefineElement(4, element, 0.0183836);
526 mix = new TGeoMixture("TPC_MYLAR", 3, 1.39);
527 element = table->GetElement(1);
528 mix->DefineElement(0, element, 0.0416667);
529 element = table->GetElement(6);
530 mix->DefineElement(1, element, 0.625);
531 element = table->GetElement(8);
532 mix->DefineElement(2, element, 0.333333);
534 case 18: // SPDBUS(AL+KPT+EPOX) - unknown composition
535 mix = new TGeoMixture("ITS_SPDBUS", 1, 1.906);
536 element = table->GetElement(9);
537 mix->DefineElement(0, element, 1.);
540 case 19: // SDD/SSD rings - unknown composition
541 mix = new TGeoMixture("ITS_SDDRINGS", 1, 1.8097);
542 element = table->GetElement(6);
543 mix->DefineElement(0, element, 1.);
546 case 20: // SPD end ladder - unknown composition
547 mix = new TGeoMixture("ITS_SPDEL", 1, 3.6374);
548 element = table->GetElement(22);
549 mix->DefineElement(0, element, 1.);
552 case 21: // SDD end ladder - unknown composition
553 mix = new TGeoMixture("ITS_SDDEL", 1, 0.3824);
554 element = table->GetElement(30);
555 mix->DefineElement(0, element, 1.);
558 case 22: // SSD end ladder - unknown composition
559 mix = new TGeoMixture("ITS_SSDEL", 1, 0.68);
560 element = table->GetElement(16);
561 mix->DefineElement(0, element, 1.);
566 printf("Patched with mixture %s\n", mix->GetName());
570 //_____________________________________________________________________________
571 void TFlukaMCGeometry::CreateFlukaMatFile(const char *fname)
573 // ==== from FLUGG ====
574 // NAMES OF ELEMENTS AND COMPOUNDS: the names must be written in upper case,
575 // according to the fluka standard. In addition,. they must be equal to the
576 // names of the fluka materials - see fluka manual - in order that the
577 // program load the right cross sections, and equal to the names included in
578 // the .pemf. Otherwise the user must define the LOW-MAT CARDS, and make his
579 // own .pemf, in order to get the right cross sections loaded in memory.
580 // Materials defined by FLUKA
582 gGeoManager->Export("geometry.root");
583 if (fname) sname = fname;
584 else sname = "flukaMat.inp";
586 out.open(sname.Data(), ios::out);
588 Fatal("CreateFlukaMatFile", "could not open file %s for writing", sname.Data());
591 PrintHeader(out, "MATERIALS AND COMPOUNDS");
592 PrintHeader(out, "MATERIALS");
594 Int_t counttothree, nelem;
596 TGeoElementTable *table = gGeoManager->GetElementTable();
597 TGeoElement *element;
598 element = table->GetElement(13);
599 element->SetTitle("ALUMINUM"); // this is how FLUKA likes it ...
600 element = table->GetElement(15);
601 element->SetTitle("PHOSPHO"); // same story ...
602 Int_t nelements = table->GetNelements();
603 TList *matlist = gGeoManager->GetListOfMaterials();
605 Int_t nmater = matlist->GetSize();
607 Int_t newind = 26; // here non predefined materials start
610 TGeoMixture *mix = 0;
613 // Create all needed elements
614 for (Int_t i=1; i<nelements; i++) {
615 element = table->GetElement(i);
616 // skip elements which are not defined
617 if (!element->IsUsed() && !element->IsDefined()) continue;
618 matname = element->GetTitle();
619 ToFlukaString(matname);
622 mat = new TGeoMaterial(matname, element->A(), element->Z(), rho);
623 // Check if the element has been predefined by FLUKA
624 Int_t pmid = GetPredefinedMaterialId(Int_t(element->Z()));
628 mat->SetIndex(newind++);
633 objstr = new TObjString(matname.Data());
634 fMatNames->Add(objstr);
639 // Adjust material names and add them to FLUKA list
640 for (i=0; i<nmater; i++) {
641 mat = (TGeoMaterial*)matlist->At(i);
642 if (!mat->IsUsed()) continue;
645 rho = mat->GetDensity();
646 if (mat->GetZ()<0.001) {
647 mat->SetIndex(2); // vacuum, built-in inside FLUKA
650 matname = mat->GetName();
651 FlukaMatName(matname);
653 mat->SetIndex(newind++);
654 objstr = new TObjString(matname.Data());
656 fMatNames->Add(objstr);
660 // Dump all elements with MATERIAL cards
661 for (i=0; i<nfmater; i++) {
662 mat = (TGeoMaterial*)fMatList->At(i);
663 // mat->SetUsed(kFALSE);
665 out << setw(10) << "MATERIAL ";
666 out.setf(static_cast<std::ios::fmtflags>(0),std::ios::floatfield);
667 objstr = (TObjString*)fMatNames->At(i);
668 matname = objstr->GetString();
671 rho = mat->GetDensity();
672 if (mat->IsMixture()) {
673 out << setw(10) << " ";
674 out << setw(10) << " ";
675 mix = (TGeoMixture*)mat;
677 out << setw(10) << setiosflags(ios::fixed) << setprecision(1) << z;
678 out << setw(10) << setprecision(3) << a;
680 out.setf(static_cast<std::ios::fmtflags>(0),std::ios::floatfield);
681 out << setw(10) << setiosflags(ios::scientific) << setprecision(3) << rho;
682 out.setf(static_cast<std::ios::fmtflags>(0),std::ios::floatfield);
683 out << setw(10) << setiosflags(ios::fixed) << setprecision(1) << Double_t(mat->GetIndex());
684 out << setw(10) << " ";
685 out << setw(10) << " ";
686 out << setw(8) << matname.Data() << endl;
688 // add LOW-MAT card for NEON to associate with ARGON neutron xsec
690 out << setw(10) << "LOW-MAT ";
691 out.setf(static_cast<std::ios::fmtflags>(0),std::ios::floatfield);
692 out << setw(10) << setiosflags(ios::fixed) << setprecision(1) << Double_t(mat->GetIndex());
693 out << setw(10) << setiosflags(ios::fixed) << setprecision(1) << 18.;
694 out << setw(10) << setiosflags(ios::fixed) << setprecision(1) << -2.;
695 out << setw(10) << setiosflags(ios::fixed) << setprecision(1) << 293.;
696 out << setw(10) << " ";
697 out << setw(10) << " ";
698 // out << setw(8) << matname.Data() << endl;
699 out << setw(8) << " " << endl;
702 element = table->GetElement((int)z);
703 TString elename = element->GetTitle();
704 ToFlukaString(elename);
705 if( matname.CompareTo( elename ) != 0 ) {
706 out << setw(10) << "LOW-MAT ";
707 out.setf(static_cast<std::ios::fmtflags>(0),std::ios::floatfield);
708 out << setw(10) << setiosflags(ios::fixed) << setprecision(1) << Double_t(mat->GetIndex());
709 out << setw(10) << setiosflags(ios::fixed) << setprecision(1) << z;
710 out << setw(10) << setiosflags(ios::fixed) << setprecision(1) << " ";
711 out << setw(10) << setiosflags(ios::fixed) << setprecision(1) << " ";
712 out << setw(10) << " ";
713 out << setw(10) << " ";
714 // missing material at Low Energy Cross Section Table
715 if( (int)z==10 || (int)z==21 || (int)z==34 || (int)z==37 || (int)z==39 || (int)z==44 ||
716 (int)z==45 || (int)z==46 || (int)z==52 || (int)z==57 || (int)z==59 || (int)z==60 ||
717 (int)z==61 || (int)z==65 || (int)z==66 || (int)z==67 || (int)z==68 || (int)z==69 ||
718 (int)z==70 || (int)z==71 || (int)z==72 || (int)z==76 || (int)z==77 || (int)z==78 ||
719 (int)z==81 || (int)z==84 || (int)z==85 || (int)z==86 || (int)z==87 || (int)z==88 ||
720 (int)z==89 || (int)z==91 )
721 out << setw(8) << "UNKNOWN " << endl;
723 out << setw(8) << elename.Data() << endl;
724 // out << setw(8) << " " << endl;
730 out << setw(10) << "COMPOUND ";
731 nelem = mix->GetNelements();
732 objstr = (TObjString*)fMatNames->At(i);
733 matname = objstr->GetString();
734 for (j=0; j<nelem; j++) {
735 w = (mix->GetWmixt())[j];
736 if (w<0.00001) w=0.00001;
737 z = (mix->GetZmixt())[j];
738 a = (mix->GetAmixt())[j];
739 idmat = GetElementIndex(Int_t(z));
740 if (!idmat) Error("CreateFlukaMatFile", "element with Z=%f not found", z);
741 out.setf(static_cast<std::ios::fmtflags>(0),std::ios::floatfield);
742 out << setw(10) << setiosflags(ios::fixed) << setprecision(6) << -w;
743 out.setf(static_cast<std::ios::fmtflags>(0),std::ios::floatfield);
744 out << setw(10) << setiosflags(ios::fixed) << setprecision(1) << Double_t(idmat);
746 if (counttothree == 3) {
747 out << matname.Data();
749 if ( (j+1) != nelem) out << setw(10) << "COMPOUND ";
754 for (j=0; j<(3-(nelem%3)); j++)
755 out << setw(10) << " " << setw(10) << " ";
756 out << matname.Data();
760 Int_t nvols = gGeoManager->GetListOfUVolumes()->GetEntriesFast()-1;
762 // Now print the material assignments
763 Double_t flagfield = 0.;
764 printf("#############################################################\n");
765 if (gFluka->IsFieldEnabled()) {
767 printf("Magnetic field enabled\n");
768 } else printf("Magnetic field disabled\n");
769 printf("#############################################################\n");
771 PrintHeader(out, "TGEO MATERIAL ASSIGNMENTS");
772 for (i=1; i<=nvols; i++) {
774 vol = gGeoManager->GetVolume(i);
775 if ((med = vol->GetMedium()) == 0) continue;
776 mat = med->GetMaterial();
777 idmat = mat->GetIndex();
778 for (Int_t j=0; j<nfmater; j++) {
779 mat = (TGeoMaterial*)fMatList->At(j);
780 if (mat->GetIndex() == idmat) mat->SetUsed(kTRUE);
783 Float_t hasfield = (vol->GetMedium()->GetParam(1) > 0) ? flagfield : 0.;
784 out << "* Assigning material: " << vol->GetMedium()->GetMaterial()->GetName() << " to Volume: " << vol->GetName();
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";
792 out << setw(10) << "0.0";
793 out << setw(10) << setiosflags(ios::fixed) << hasfield;
794 out << setw(10) << "0.0";
799 fDummyRegion = nvols+1;
800 out << "* Dummy region: " << endl;
801 out << setw(10) << "ASSIGNMAT ";
802 out.setf(static_cast<std::ios::fmtflags>(0),std::ios::floatfield);
803 out << setw(10) << setiosflags(ios::fixed) << idmat;
804 out << setw(10) << setiosflags(ios::fixed) << fDummyRegion;
805 out << setw(10) << "0.0";
806 out << setw(10) << "0.0";
807 out << setw(10) << "0.0";
808 out << setw(10) << "0.0" << endl;
810 // fLastMaterial = nfmater+2;
811 fLastMaterial = newind;
814 //_____________________________________________________________________________
815 void TFlukaMCGeometry::PrintHeader(ofstream &out, const char *text) const
817 // Print a FLUKA header.
818 out << "*\n" << "*\n" << "*\n";
819 out << "********************* " << text << " *********************\n"
821 out << "*...+....1....+....2....+....3....+....4....+....5....+....6....+....7..."
826 //_____________________________________________________________________________
827 Int_t TFlukaMCGeometry::RegionId() const
829 // Returns current region id <-> TGeo node id
830 if (gGeoManager->IsOutside()) return 0;
831 return gGeoManager->GetCurrentNode()->GetUniqueID();
834 //_____________________________________________________________________________
835 Int_t TFlukaMCGeometry::GetElementIndex(Int_t z) const
837 // Get index of a material having a given Z element.
838 TIter next(fMatList);
841 while ((mat=(TGeoMaterial*)next())) {
842 if (mat->IsMixture()) continue;
843 if (mat->GetElement()->Z() == z) return mat->GetIndex();
848 //_____________________________________________________________________________
849 void TFlukaMCGeometry::SetMreg(Int_t mreg, Int_t lttc)
851 // Update if needed next history;
852 // if (gFluka->GetDummyBoundary()==2) {
853 // gGeoManager->CdNode(fNextLattice-1);
856 if (lttc == TFlukaMCGeometry::kLttcOutside) {
857 fCurrentRegion = NofVolumes()+2;
858 fCurrentLattice = lttc;
859 gGeoManager->CdTop();
860 gGeoManager->SetOutside(kTRUE);
862 if (lttc == TFlukaMCGeometry::kLttcVirtual) return;
864 Error("TFlukaMCGeometry::SetMreg","Invalide lattice %i",lttc);
867 fCurrentRegion = mreg;
868 fCurrentLattice = lttc;
870 Int_t crtlttc = gGeoManager->GetCurrentNodeId()+1;
871 if (crtlttc == lttc) return;
872 gGeoManager->CdNode(lttc-1);
873 while (gGeoManager->GetCurrentVolume()->IsAssembly()) gGeoManager->CdUp();
876 //_____________________________________________________________________________
877 void TFlukaMCGeometry::SetCurrentRegion(Int_t mreg, Int_t latt)
879 // Set index/history for next entered region
880 fCurrentRegion = mreg;
881 fCurrentLattice = latt;
884 //_____________________________________________________________________________
885 void TFlukaMCGeometry::SetNextRegion(Int_t mreg, Int_t latt)
887 // Set index/history for next entered region
892 //_____________________________________________________________________________
893 void TFlukaMCGeometry::ToFlukaString(TString &str) const
895 // ToFlukaString converts an string to something usefull in FLUKA:
898 // * Replace ' ' by '_'
899 if (str.Length()<8) {
904 for (ilast=7; ilast>0; ilast--) if (str(ilast)!=' ') break;
906 for (Int_t pos=0; pos<ilast; pos++)
907 if (str(pos)==' ') str.Replace(pos,1,"_",1);
911 //_____________________________________________________________________________
912 void TFlukaMCGeometry::FlukaMatName(TString &str) const
914 // Strip the detector name
915 TObjArray * tokens = str.Tokenize("_");
916 Int_t ntok = tokens->GetEntries();
918 TString head = ((TObjString*) tokens->At(0))->GetString();
919 Int_t nhead = head.Length();
920 str = str.Remove(0, nhead + 1);
924 // Convert a name to upper case 8 chars.
927 for (ilast=7; ilast>0; ilast--) if (str(ilast)!=' ') break;
928 if (ilast>5) ilast = 5;
930 TIter next(fMatNames);
934 while ((objstr=(TObjString*)next())) {
935 matname = objstr->GetString();
936 if (matname == str) {
940 sprintf(&number[1], "%d", index);
941 } else if (index<100) {
942 sprintf(number, "%d", index);
944 Error("FlukaMatName", "Too many materials %s", str.Data());
947 str.Replace(ilast+1, 2, number);
953 //______________________________________________________________________________
954 void TFlukaMCGeometry::Vname(const char *name, char *vname) const
957 // convert name to upper case. Make vname at least 4 chars
959 Int_t l = strlen(name);
962 for (i=0;i<l;i++) vname[i] = toupper(name[i]);
963 for (i=l;i<4;i++) vname[i] = ' ';
968 //______________________________________________________________________________
969 Int_t TFlukaMCGeometry::GetNstep()
971 // return gNstep for debug propose
976 Int_t TFlukaMCGeometry::GetPredefinedMaterialId(Int_t z) const
978 // Get predifined material id from Z if present in list
979 const Int_t kMax = 25;
981 static Int_t idFluka[kMax] =
982 {-1, -1, 1, 2, 4, 6, 7, 8, 12, 13,
983 26, 29, 47, 14, 79, 80, 82, 73, 11, 18,
988 for (Int_t i = 0; i < kMax; i++)
990 if (z == idFluka[i]) {
1000 // FLUKA GEOMETRY WRAPPERS - to replace FLUGG wrappers
1002 //_____________________________________________________________________________
1003 Int_t idnrwr(const Int_t & /*nreg*/, const Int_t & /*mlat*/)
1006 // Wrapper for setting DNEAR option on fluka side. Must return 0
1007 // if user doesn't want Fluka to use DNEAR to compute the
1008 // step (the same effect is obtained with the GLOBAL (WHAT(3)=-1)
1009 // card in fluka input), returns 1 if user wants Fluka always to
1010 // use DNEAR (in this case, be sure that GEANT4 DNEAR is unique,
1011 // coming from all directions!!!)
1012 if (gMCGeom->IsDebugging()) printf("========== Dummy IDNRWR\n");
1016 //_____________________________________________________________________________
1017 void g1wr(Double_t &pSx, Double_t &pSy, Double_t &pSz,
1018 Double_t *pV, Int_t &oldReg , const Int_t &oldLttc, Double_t &propStep,
1019 Int_t &/*nascFlag*/, Double_t &retStep, Int_t &newReg,
1020 Double_t &saf, Int_t &newLttc, Int_t <tcFlag,
1021 Double_t *sLt, Int_t *jrLt)
1024 // Initialize FLUKa point and direction;
1025 static Int_t ierr = 0;
1027 // gMCGeom->SetDebugMode(kTRUE);
1033 Int_t olttc = oldLttc;
1035 gGeoManager->FindNode(pSx,pSy,pSz);
1036 olttc = gGeoManager->GetCurrentNodeId()+1;
1037 oldReg = gGeoManager->GetCurrentVolume()->GetNumber();
1040 if (gMCGeom->IsDebugging()) {
1041 cout << "g1wr gNstep=" << gNstep << " oldReg="<< oldReg <<" olttc="<< olttc
1042 << " track=" << TRACKR.ispusr[mkbmx2-1] << endl;
1043 cout << " point: (" << pSx << ", " << pSy << ", " << pSz << ") dir: ("
1044 << pV[0] << ", " << pV[1] << ", " << pV[2] << ")" << endl;
1048 Int_t ccreg=0,cclat=0;
1049 gMCGeom->GetCurrentRegion(ccreg,cclat);
1050 Bool_t crossed = (ccreg==oldReg && cclat==olttc)?kFALSE:kTRUE;
1052 gMCGeom->SetCurrentRegion(oldReg, olttc);
1053 // Initialize default return values
1055 jrLt[lttcFlag] = olttc;
1056 sLt[lttcFlag] = propStep;
1057 jrLt[lttcFlag+1] = -1;
1058 sLt[lttcFlag+1] = 0.;
1061 Bool_t crossedDummy = (oldReg == gFluka->GetDummyRegion())?kTRUE:kFALSE;
1062 Int_t curLttc, curReg;
1064 // FLUKA crossed the dummy boundary - update new region/history
1065 retStep = TGeoShape::Tolerance();
1067 gMCGeom->GetNextRegion(newReg, newLttc);
1068 gMCGeom->SetMreg(newReg, newLttc);
1069 sLt[lttcFlag] = TGeoShape::Tolerance(); // null step in current region
1071 jrLt[lttcFlag] = newLttc;
1072 sLt[lttcFlag] = TGeoShape::Tolerance(); // null step in next region
1073 jrLt[lttcFlag+1] = -1;
1074 sLt[lttcFlag+1] = 0.; // null step in next region;
1075 if (gMCGeom->IsDebugging()) printf("=> crossed dummy!! newReg=%i newLttc=%i\n", newReg, newLttc);
1079 // Reset outside flag
1080 gGeoManager->SetOutside(kFALSE);
1082 curLttc = gGeoManager->GetCurrentNodeId()+1;
1083 curReg = gGeoManager->GetCurrentVolume()->GetNumber();
1084 if (olttc != curLttc) {
1085 // FLUKA crossed the boundary : we trust that the given point is really there,
1086 // so we just update TGeo state
1087 gGeoManager->CdNode(olttc-1);
1088 curLttc = gGeoManager->GetCurrentNodeId()+1;
1089 curReg = gGeoManager->GetCurrentVolume()->GetNumber();
1091 // Now the current TGeo state reflects the FLUKA state
1093 gGeoManager->SetCurrentPoint(pSx, pSy, pSz);
1094 gGeoManager->SetCurrentDirection(pV);
1095 Double_t pt[3], local[3], ldir[3];
1096 Int_t pid = TRACKR.jtrack;
1100 gGeoManager->MasterToLocal(pt,local);
1101 gGeoManager->MasterToLocalVect(pV,ldir);
1103 Bool_t valid = gGeoManager->GetCurrentVolume()->Contains(local);
1105 printf("location not valid in %s pid=%i\n", gGeoManager->GetPath(),pid);
1106 printf("local:(%f, %f, %f) ldir:(%f, %f, %f)\n", local[0],local[1],local[2],ldir[0],ldir[1],ldir[2]);
1107 // gGeoManager->FindNode();
1108 // printf(" -> actual location: %s\n", gGeoManager->GetPath());
1111 Double_t pstep = propStep;
1112 Double_t snext = propStep;
1113 const Double_t epsil = 0.9999999 * TGeoShape::Tolerance();
1114 // This should never happen !!!
1115 if (pstep<TGeoShape::Tolerance()) {
1116 printf("Proposed step is 0 !!!\n");
1117 pstep = 2.*TGeoShape::Tolerance();
1120 gGeoManager->FindNextBoundaryAndStep(pstep);
1121 snext = gGeoManager->GetStep();
1123 if (snext <= TGeoShape::Tolerance()) {
1124 // printf("FLUKA crossed boundary but snext=%g\n", snext);
1128 snext += TGeoShape::Tolerance();
1132 gGeoManager->FindNextBoundaryAndStep(pstep, kTRUE);
1133 snext = gGeoManager->GetStep();
1134 saf = gGeoManager->GetSafeDistance();
1135 if (snext <= TGeoShape::Tolerance()) {
1136 // printf("FLUKA put particle on bondary without crossing\n");
1141 snext += TGeoShape::Tolerance();
1148 // printf("%d snext=%g\n", ierr, snext);
1151 // printf("Too many null steps - sending error code -33...\n");
1152 newReg = -33; // Error code
1158 NORLAT.distn = snext;
1159 NORLAT.xn[0] += snext*pV[0];
1160 NORLAT.xn[1] += snext*pV[1];
1161 NORLAT.xn[2] += snext*pV[2];
1162 if (!gGeoManager->IsOnBoundary()) {
1163 // Next boundary further than proposed step, which is approved
1164 if (saf>propStep) saf = propStep;
1166 sLt[lttcFlag] = propStep;
1169 if (saf>snext) saf = snext; // Safety should be less than the proposed step if a boundary will be crossed
1170 gGeoManager->SetCurrentPoint(pSx,pSy,pSz);
1171 newLttc = (gGeoManager->IsOutside())?(TFlukaMCGeometry::kLttcOutside):gGeoManager->GetCurrentNodeId()+1;
1172 newReg = (gGeoManager->IsOutside())?(gMCGeom->NofVolumes()+2):gGeoManager->GetCurrentVolume()->GetNumber();
1173 if (gMCGeom->IsDebugging()) printf("=> newReg=%i newLttc=%i\n", newReg, newLttc);
1175 // We really crossed the boundary, but is it the same region ?
1176 gMCGeom->SetNextRegion(newReg, newLttc);
1178 if ( ((newReg==oldReg && newLttc!=olttc) || (oldReg!=newReg && olttc==newLttc) ) && pid!=-1) {
1179 // Virtual boundary between replicants
1180 newReg = gFluka->GetDummyRegion();
1181 newLttc = TFlukaMCGeometry::kLttcVirtual;
1182 if (gMCGeom->IsDebugging()) printf("=> virtual boundary!! newReg=%i newLttc=%i\n", newReg, newLttc);
1186 sLt[lttcFlag] = snext;
1188 jrLt[lttcFlag] = newLttc;
1189 sLt[lttcFlag] = snext;
1190 jrLt[lttcFlag+1] = -1;
1192 sLt[lttcFlag+1] = 0.;
1193 gGeoManager->SetOutside(kFALSE);
1194 gGeoManager->CdNode(olttc-1);
1195 if (gMCGeom->IsDebugging()) {
1196 printf("=> snext=%g safe=%g\n", snext, saf);
1197 for (Int_t i=0; i<lttcFlag+1; i++) printf(" jrLt[%i]=%i sLt[%i]=%g\n", i,jrLt[i],i,sLt[i]);
1201 //_____________________________________________________________________________
1205 if (gMCGeom->IsDebugging()) printf("========== Dummy G1RTWR\n");
1208 //_____________________________________________________________________________
1209 void conhwr(Int_t & intHist, Int_t & incrCount)
1211 if (gMCGeom->IsDebugging()) printf("========== Dummy CONHWR intHist=%d incrCount=%d currentNodeId=%d\n",
1212 intHist, incrCount, gGeoManager->GetCurrentNodeId()+1 );
1213 // if( incrCount != -1 ) {
1214 // if (intHist==0) gGeoManager->CdTop();
1215 // else gGeoManager->CdNode(intHist-1);
1217 // intHist = gGeoManager->GetCurrentNodeId()+1;
1220 //_____________________________________________________________________________
1221 void inihwr(Int_t &intHist)
1223 if (gMCGeom->IsDebugging())
1224 printf("========== Inside INIHWR -> reinitializing history: %i \n", intHist);
1225 if (gGeoManager->IsOutside()) gGeoManager->CdTop();
1227 // printf("=== wrong history number\n");
1230 if (intHist==0) gGeoManager->CdTop();
1231 else gGeoManager->CdNode(intHist-1);
1232 if (gMCGeom->IsDebugging()) {
1233 printf(" --- current path: %s\n", gGeoManager->GetPath());
1234 printf("<= INIHWR\n");
1238 //_____________________________________________________________________________
1239 void jomiwr(const Int_t & /*nge*/, const Int_t & /*lin*/, const Int_t & /*lou*/,
1242 // Geometry initialization wrapper called by FLUKAM. Provides to FLUKA the
1243 // number of regions (volumes in TGeo)
1244 // build application geometry
1245 if (gMCGeom->IsDebugging()) printf("========== Inside JOMIWR\n");
1246 flukaReg = gGeoManager->GetListOfUVolumes()->GetEntriesFast()+1;
1247 if (gMCGeom->IsDebugging()) printf("<= JOMIWR: last region=%i\n", flukaReg);
1250 //_____________________________________________________________________________
1251 void lkdbwr(Double_t &pSx, Double_t &pSy, Double_t &pSz,
1252 Double_t *pV, const Int_t &oldReg, const Int_t &oldLttc,
1253 Int_t &flagErr, Int_t &newReg, Int_t &newLttc)
1255 if (gMCGeom->IsDebugging()) {
1256 printf("========== Inside LKDBWR (%f, %f, %f)\n",pSx, pSy, pSz);
1257 printf(" in: pV=(%f, %f, %f)\n", pV[0], pV[1], pV[2]);
1258 printf(" in: oldReg=%i oldLttc=%i\n", oldReg, oldLttc);
1260 lkwr(pSx,pSy,pSz,pV,oldReg,oldLttc,flagErr,newReg,newLttc);
1263 //_____________________________________________________________________________
1264 void lkfxwr(Double_t &pSx, Double_t &pSy, Double_t &pSz,
1265 Double_t *pV, const Int_t &oldReg, const Int_t &oldLttc,
1266 Int_t &flagErr, Int_t &newReg, Int_t &newLttc)
1268 if (gMCGeom->IsDebugging()) {
1269 printf("========== Inside LKFXWR (%f, %f, %f)\n",pSx, pSy, pSz);
1270 printf(" in: pV=(%f, %f, %f)\n", pV[0], pV[1], pV[2]);
1271 printf(" in: oldReg=%i oldLttc=%i\n", oldReg, oldLttc);
1273 lkwr(pSx,pSy,pSz,pV,oldReg,oldLttc,flagErr,newReg,newLttc);
1276 //_____________________________________________________________________________
1277 void lkmgwr(Double_t &pSx, Double_t &pSy, Double_t &pSz,
1278 Double_t *pV, const Int_t &oldReg, const Int_t &oldLttc,
1279 Int_t &flagErr, Int_t &newReg, Int_t &newLttc)
1281 if (gMCGeom->IsDebugging()) {
1282 printf("========== Inside LKMGWR (%f, %f, %f)\n",pSx, pSy, pSz);
1283 printf(" in: pV=(%f, %f, %f)\n", pV[0], pV[1], pV[2]);
1284 printf(" in: oldReg=%i oldLttc=%i\n", oldReg, oldLttc);
1286 lkwr(pSx,pSy,pSz,pV,oldReg,oldLttc,flagErr,newReg,newLttc);
1289 //_____________________________________________________________________________
1290 void lkwr(Double_t &pSx, Double_t &pSy, Double_t &pSz,
1291 Double_t *pV, const Int_t &oldReg, const Int_t &oldLttc,
1292 Int_t &flagErr, Int_t &newReg, Int_t &newLttc)
1294 if (gMCGeom->IsDebugging()) {
1295 printf("========== Inside LKWR (%f, %f, %f)\n",pSx, pSy, pSz);
1296 printf(" in: pV=(%f, %f, %f)\n", pV[0], pV[1], pV[2]);
1297 printf(" in: oldReg=%i oldLttc=%i\n", oldReg, oldLttc);
1300 Double_t epsil = 1000.*TGeoShape::Tolerance();
1301 TGeoNode *node = gGeoManager->FindNode(pSx+epsil*pV[0], pSy+epsil*pV[1], pSz+epsil*pV[2]);
1302 if (gGeoManager->IsOutside()) {
1303 newReg = gMCGeom->NofVolumes()+2;
1304 newLttc = TFlukaMCGeometry::kLttcOutside;
1305 gGeoManager->SetOutside(kFALSE);
1306 if (oldLttc>0 && oldLttc<newLttc) gGeoManager->CdNode(oldLttc-1);
1309 gGeoManager->SetOutside(kFALSE);
1310 newReg = node->GetVolume()->GetNumber();
1311 newLttc = gGeoManager->GetCurrentNodeId()+1;
1312 if (oldLttc==TFlukaMCGeometry::kLttcOutside || oldLttc==0) return;
1314 Int_t dummy = gFluka->GetDummyRegion();
1315 if (oldReg==dummy) {
1316 Int_t newreg1, newlttc1;
1317 gMCGeom->GetNextRegion(newreg1, newlttc1);
1318 if (newreg1==newReg && newlttc1==newLttc) {
1320 newLttc = TFlukaMCGeometry::kLttcVirtual;
1322 if (gMCGeom->IsDebugging()) printf(" virtual boundary (oldReg==dummy) !! newReg=%i newLttc=%i\n", newReg, newLttc);
1327 if (oldReg==newReg && oldLttc!=newLttc) {
1329 newLttc = TFlukaMCGeometry::kLttcVirtual;
1330 if (gMCGeom->IsDebugging()) printf(" virtual boundary!! newReg=%i newLttc=%i\n", newReg, newLttc);
1333 if( oldReg!=newReg && oldLttc==newLttc ) {
1334 // this should not happen!! ??? Ernesto
1335 // cout << " lkwr oldReg!=newReg ("<< oldReg <<"!="<< newReg
1336 // << ") && oldLttc==newLttc ("<< newLttc <<") !!!!!!!!!!!!!!!!!" << endl;
1338 newLttc = TFlukaMCGeometry::kLttcVirtual;
1342 if (gMCGeom->IsDebugging()) {
1343 printf(" LKWR: newReg=%i newLttc=%i\n", newReg, newLttc);
1347 //_____________________________________________________________________________
1348 void nrmlwr(Double_t &pSx, Double_t &pSy, Double_t &pSz,
1349 Double_t &pVx, Double_t &pVy, Double_t &pVz,
1350 Double_t *norml, const Int_t &oldReg,
1351 const Int_t &newReg, Int_t &flagErr)
1353 if (gMCGeom->IsDebugging()) {
1354 printf("========== Inside NRMLWR (%g, %g, %g, %g, %g, %g)\n", pSx,pSy,pSz,pVx,pVy,pVz);
1355 printf(" (%g, %g, %g)\n", NORLAT.xn[0], NORLAT.xn[1], NORLAT.xn[2]);
1356 printf(" oldReg=%i, newReg=%i\n", oldReg,newReg);
1358 gGeoManager->SetCurrentPoint(NORLAT.xn[0], NORLAT.xn[1], NORLAT.xn[2]);
1359 gGeoManager->SetCurrentDirection(pVx, pVy, pVz);
1360 Double_t *dnorm = gGeoManager->FindNormalFast();
1363 printf(" ERROR: Cannot compute fast normal\n");
1369 norml[0] = -dnorm[0];
1370 norml[1] = -dnorm[1];
1371 norml[2] = -dnorm[2];
1374 if (gMCGeom->IsDebugging()) {
1375 printf(" normal to boundary: (%g, %g, %g)\n", norml[0], norml[1], norml[2]);
1376 printf("<= NRMLWR\n");
1381 //_____________________________________________________________________________
1382 void rgrpwr(const Int_t & /*flukaReg*/, const Int_t & /*ptrLttc*/, Int_t & /*g4Reg*/,
1383 Int_t * /*indMother*/, Int_t * /*repMother*/, Int_t & /*depthFluka*/)
1385 if (gMCGeom->IsDebugging()) printf("=> Dummy RGRPWR\n");
1388 //_____________________________________________________________________________
1389 Int_t isvhwr(const Int_t &check, const Int_t & intHist)
1392 // Wrapper for saving current navigation history (fCheck=default)
1393 // and returning its pointer. If fCheck=-1 copy of history pointed
1394 // by intHist is made in NavHistWithCount object, and its pointer
1395 // is returned. fCheck=1 and fCheck=2 cases are only in debugging
1396 // version: an array is created by means of FGeometryInit functions
1397 // (but could be a static int * ptrArray = new int[10000] with
1398 // file scope as well) that stores a flag for deleted/undeleted
1399 // histories and at the end of event is checked to verify that
1400 // all saved history objects have been deleted.
1402 // For TGeo, just return the current node ID. No copy need to be made.
1404 if (gMCGeom->IsDebugging()) printf("=> Inside ISVHWR check=%d intHist=%d\n", check, intHist);
1405 if (check<0) return intHist;
1406 Int_t histInt = gGeoManager->GetCurrentNodeId()+1;
1407 if (gMCGeom->IsDebugging()) printf("<= ISVHWR: history is: %i in: %s\n", histInt, gGeoManager->GetPath());