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"
29 #include "TFlukaMCGeometry.h"
30 #include "TFlukaConfigOption.h"
31 #include "TGeoManager.h"
32 #include "TGeoVolume.h"
33 #include "TObjString.h"
38 # define idnrwr idnrwr_
40 # define g1rtwr g1rtwr_
41 # define conhwr conhwr_
42 # define inihwr inihwr_
43 # define jomiwr jomiwr_
44 # define lkdbwr lkdbwr_
45 # define lkfxwr lkfxwr_
46 # define lkmgwr lkmgwr_
48 # define magfld magfld_
49 # define nrmlwr nrmlwr_
50 # define rgrpwr rgrpwr_
51 # define isvhwr isvhwr_
55 # define idnrwr IDNRWR
57 # define g1rtwr G1RTWR
58 # define conhwr CONHWR
59 # define inihwr INIHWR
60 # define jomiwr JOMIWR
61 # define lkdbwr LKDBWR
62 # define lkfxwr LKFXWR
63 # define lkmgwr LKMGWR
65 # define magfld MAGFLD
66 # define nrmlwr NRMLWR
67 # define rgrpwr RGRPWR
68 # define isvhwr ISVHWR
72 //____________________________________________________________________________
76 // Prototypes for FLUKA navigation methods
78 Int_t type_of_call idnrwr(const Int_t & /*nreg*/, const Int_t & /*mlat*/);
79 void type_of_call g1wr(Double_t & /*pSx*/, Double_t & /*pSy*/, Double_t & /*pSz*/,
80 Double_t * /*pV*/, Int_t & /*oldReg*/ , const Int_t & /*oldLttc*/, Double_t & /*propStep*/,
81 Int_t & /*nascFlag*/, Double_t & /*retStep*/, Int_t & /*newReg*/,
82 Double_t & /*saf*/, Int_t & /*newLttc*/, Int_t & /*LttcFlag*/,
83 Double_t *s /*Lt*/, Int_t * /*jrLt*/);
85 void type_of_call g1rtwr();
86 void type_of_call conhwr(Int_t & /*intHist*/, Int_t * /*incrCount*/);
87 void type_of_call inihwr(Int_t & /*intHist*/);
88 void type_of_call jomiwr(const Int_t & /*nge*/, const Int_t & /*lin*/, const Int_t & /*lou*/,
89 Int_t & /*flukaReg*/);
90 void type_of_call lkdbwr(Double_t & /*pSx*/, Double_t & /*pSy*/, Double_t & /*pSz*/,
91 Double_t * /*pV*/, const Int_t & /*oldReg*/, const Int_t & /*oldLttc*/,
92 Int_t & /*flagErr*/, Int_t & /*newReg*/, Int_t & /*newLttc*/);
93 void type_of_call lkfxwr(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 lkmgwr(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 lkwr(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 magfld(const Double_t & /*pX*/, const Double_t & /*pY*/, const Double_t & /*pZ*/,
103 // Double_t & /*cosBx*/, Double_t & /*cosBy*/, Double_t & /*cosBz*/,
104 // Double_t & /*Bmag*/, Int_t & /*reg*/, Int_t & /*idiscflag*/);
105 void type_of_call nrmlwr(Double_t & /*pSx*/, Double_t & /*pSy*/, Double_t & /*pSz*/,
106 Double_t & /*pVx*/, Double_t & /*pVy*/, Double_t & /*pVz*/,
107 Double_t * /*norml*/, const Int_t & /*oldReg*/,
108 const Int_t & /*newReg*/, Int_t & /*flagErr*/);
109 void type_of_call rgrpwr(const Int_t & /*flukaReg*/, const Int_t & /*ptrLttc*/, Int_t & /*g4Reg*/,
110 Int_t * /*indMother*/, Int_t * /*repMother*/, Int_t & /*depthFluka*/);
111 Int_t type_of_call isvhwr(const Int_t & /*fCheck*/, const Int_t & /*intHist*/);
114 // TFluka global pointer
116 TFlukaMCGeometry *gMCGeom = 0;
119 ClassImp(TFlukaMCGeometry)
121 TFlukaMCGeometry* TFlukaMCGeometry::fgInstance= NULL;
123 //_____________________________________________________________________________
124 TFlukaMCGeometry::TFlukaMCGeometry(const char *name, const char *title)
125 : TNamed(name, title)
128 // Standard constructor
139 gFluka = (TFluka*)gMC;
142 fMatList = new TObjArray(256);
143 fMatNames = new TObjArray(256);
146 //_____________________________________________________________________________
147 TFlukaMCGeometry::TFlukaMCGeometry()
150 // Default constructor
161 gFluka = (TFluka*)gMC;
168 //_____________________________________________________________________________
169 TFlukaMCGeometry::~TFlukaMCGeometry()
175 if (fRegionList) delete [] fRegionList;
176 if (fMatList) delete fMatList;
177 if (fMatNames) {fMatNames->Delete(); delete fMatNames;}
178 if (gGeoManager) delete gGeoManager;
184 //_____________________________________________________________________________
185 TFlukaMCGeometry::TFlukaMCGeometry(const TFlukaMCGeometry &)
193 //_____________________________________________________________________________
194 Double_t* TFlukaMCGeometry::CreateDoubleArray(Float_t* array, Int_t size) const
196 // Converts Float_t* array to Double_t*,
197 // !! The new array has to be deleted by user.
200 Double_t* doubleArray;
202 doubleArray = new Double_t[size];
203 for (Int_t i=0; i<size; i++) doubleArray[i] = array[i];
207 doubleArray = new Double_t[1];
215 //_____________________________________________________________________________
216 Int_t TFlukaMCGeometry::GetMedium() const
218 // Get current medium number
220 TGeoNode *node = gGeoManager->GetCurrentNode();
221 if (!node) imed = gGeoManager->GetTopNode()->GetVolume()->GetMedium()->GetId();
222 else imed = node->GetVolume()->GetMedium()->GetId();
223 if (fDebug) printf("GetMedium=%i\n", imed);
227 //_____________________________________________________________________________
228 Int_t TFlukaMCGeometry::GetFlukaMaterial(Int_t imed) const
230 // Returns FLUKA material index for medium IMED
231 TGeoMedium *med = (TGeoMedium*)gGeoManager->GetListOfMedia()->At(imed-1);
233 Error("GetFlukaMaterial", "MEDIUM %i nor found", imed);
236 TGeoMaterial* mat = med->GetMaterial();
237 if (!mat->IsUsed()) return -1;
238 Int_t imatfl = med->GetMaterial()->GetIndex();
242 //_____________________________________________________________________________
243 Int_t *TFlukaMCGeometry::GetRegionList(Int_t imed, Int_t &nreg)
245 // Get an ordered list of regions matching a given medium number
247 if (!fRegionList) fRegionList = new Int_t[NofVolumes()+1];
248 TIter next(gGeoManager->GetListOfUVolumes());
251 while ((vol = (TGeoVolume*)next())) {
252 imedium = vol->GetMedium()->GetId();
253 if (imedium == imed) {
254 ireg = vol->GetNumber();
255 fRegionList[nreg++] = ireg;
261 //_____________________________________________________________________________
262 Int_t *TFlukaMCGeometry::GetMaterialList(Int_t imat, Int_t &nreg)
264 // Get an ordered list of regions matching a given medium number
266 if (!fRegionList) fRegionList = new Int_t[NofVolumes()+1];
267 TIter next(gGeoManager->GetListOfUVolumes());
269 Int_t imaterial, ireg;
270 while ((vol = (TGeoVolume*)next())) {
271 imaterial = vol->GetMedium()->GetMaterial()->GetIndex();
272 if (imaterial == imat) {
273 ireg = vol->GetNumber();
274 fRegionList[nreg++] = ireg;
280 //_____________________________________________________________________________
281 Int_t TFlukaMCGeometry::NofVolumes() const
284 // Return total number of volumes in the geometry
287 return gGeoManager->GetListOfUVolumes()->GetEntriesFast()-1;
290 //_____________________________________________________________________________
291 TGeoMaterial * TFlukaMCGeometry::GetMakeWrongMaterial(Double_t z)
293 // Try to replace a wrongly-defined material
294 static Double_t kz[23] = {7.3, 17.8184, 7.2167, 10.856, 8.875, 8.9, 7.177,
295 25.72, 6.2363, 7.1315, 47.7056, 10.6467, 7.8598, 2.10853, 10.6001, 9.1193,
296 15.3383, 4.55, 9.6502, 6.4561, 21.7963, 29.8246, 15.4021};
300 for (ind=0; ind<23; ind++) {
301 dz = TMath::Abs(z-kz[ind]);
305 printf("Cannot patch material with Z=%g\n", z);
308 TGeoMixture *mix = 0;
309 TGeoElement *element;
310 TGeoElementTable *table = gGeoManager->GetElementTable();
313 mix = new TGeoMixture("TPC_AIR", 4, 0.001205);
314 element = table->GetElement(6); // C
315 mix->DefineElement(0, element, 0.000124);
316 element = table->GetElement(7); // N
317 mix->DefineElement(1, element, 0.755267);
318 element = table->GetElement(8); // O
319 mix->DefineElement(2, element, 0.231781);
320 element = table->GetElement(18); // AR
321 mix->DefineElement(3, element, 0.012827);
323 case 1: //SDD SI CHIP
324 mix = new TGeoMixture("ITS_SDD_SI", 6, 2.4485);
325 element = table->GetElement(1);
326 mix->DefineElement(0, element, 0.004367771);
327 element = table->GetElement(6);
328 mix->DefineElement(1, element, 0.039730642);
329 element = table->GetElement(7);
330 mix->DefineElement(2, element, 0.001396798);
331 element = table->GetElement(8);
332 mix->DefineElement(3, element, 0.01169634);
333 element = table->GetElement(14);
334 mix->DefineElement(4, element, 0.844665);
335 element = table->GetElement(47);
336 mix->DefineElement(5, element, 0.09814344903);
339 mix = new TGeoMixture("ITS_WATER", 2, 1.0);
340 element = table->GetElement(1);
341 mix->DefineElement(0, element, 0.111898344);
342 element = table->GetElement(8);
343 mix->DefineElement(1, element, 0.888101656);
346 mix = new TGeoMixture("ITS_CERAMICS", 5, 3.6);
347 element = table->GetElement(8);
348 mix->DefineElement(0, element, 0.59956);
349 element = table->GetElement(13);
350 mix->DefineElement(1, element, 0.3776);
351 element = table->GetElement(14);
352 mix->DefineElement(2, element, 0.00933);
353 element = table->GetElement(24);
354 mix->DefineElement(3, element, 0.002);
355 element = table->GetElement(25);
356 mix->DefineElement(4, element, 0.0115);
359 mix = new TGeoMixture("MUON_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);
370 mix = new TGeoMixture("G10FR4", 4, 1.8);
371 element = table->GetElement(1);
372 mix->DefineElement(0, element, 0.19);
373 element = table->GetElement(6);
374 mix->DefineElement(1, element, 0.18);
375 element = table->GetElement(8);
376 mix->DefineElement(2, element, 0.35);
377 element = table->GetElement(14);
378 mix->DefineElement(3, element, 0.28);
381 mix = new TGeoMixture("ITS_KAPTON", 4, 1.3);
382 element = table->GetElement(1);
383 mix->DefineElement(0, element, 0.026363415);
384 element = table->GetElement(6);
385 mix->DefineElement(1, element, 0.6911272);
386 element = table->GetElement(7);
387 mix->DefineElement(2, element, 0.073271325);
388 element = table->GetElement(8);
389 mix->DefineElement(3, element, 0.209238060);
392 mix = new TGeoMixture("ITS_INOX", 9, 7.9);
393 element = table->GetElement(6);
394 mix->DefineElement(0, element, 0.0003);
395 element = table->GetElement(14);
396 mix->DefineElement(1, element, 0.01);
397 element = table->GetElement(15);
398 mix->DefineElement(2, element, 0.00045);
399 element = table->GetElement(16);
400 mix->DefineElement(3, element, 0.0003);
401 element = table->GetElement(24);
402 mix->DefineElement(4, element, 0.17);
403 element = table->GetElement(25);
404 mix->DefineElement(5, element, 0.02);
405 element = table->GetElement(26);
406 mix->DefineElement(6, element, 0.654);
407 element = table->GetElement(28);
408 mix->DefineElement(7, element, 0.12);
409 element = table->GetElement(42);
410 mix->DefineElement(8, element, 0.025);
413 mix = new TGeoMixture("ROHACELL", 4, 0.05);
414 element = table->GetElement(1);
415 mix->DefineElement(0, element, 0.07836617);
416 element = table->GetElement(6);
417 mix->DefineElement(1, element, 0.64648941);
418 element = table->GetElement(7);
419 mix->DefineElement(2, element, 0.08376983);
420 element = table->GetElement(8);
421 mix->DefineElement(3, element, 0.19137459);
424 mix = new TGeoMixture("ITS_SDD-C-AL", 5, 1.9837);
425 element = table->GetElement(1);
426 mix->DefineElement(0, element, 0.022632);
427 element = table->GetElement(6);
428 mix->DefineElement(1, element, 0.8176579);
429 element = table->GetElement(7);
430 mix->DefineElement(2, element, 0.0093488);
431 element = table->GetElement(8);
432 mix->DefineElement(3, element, 0.0503618);
433 element = table->GetElement(13);
434 mix->DefineElement(4, element, 0.1);
437 mix = new TGeoMixture("ITS_X7R-CAP", 7, 6.72);
438 element = table->GetElement(8);
439 mix->DefineElement(0, element, 0.085975822);
440 element = table->GetElement(22);
441 mix->DefineElement(1, element, 0.084755042);
442 element = table->GetElement(28);
443 mix->DefineElement(2, element, 0.038244751);
444 element = table->GetElement(29);
445 mix->DefineElement(3, element, 0.009471271);
446 element = table->GetElement(50);
447 mix->DefineElement(4, element, 0.321736471);
448 element = table->GetElement(56);
449 mix->DefineElement(5, element, 0.251639432);
450 element = table->GetElement(82);
451 mix->DefineElement(6, element, 0.2081768);
453 case 11: // SDD ruby sph. Al2O3
454 mix = new TGeoMixture("ITS_AL2O3", 2, 3.97);
455 element = table->GetElement(8);
456 mix->DefineElement(0, element, 0.5293);
457 element = table->GetElement(13);
458 mix->DefineElement(1, element, 0.4707);
460 case 12: // SDD HV microcable
461 mix = new TGeoMixture("ITS_HV-CABLE", 5, 1.6087);
462 element = table->GetElement(1);
463 mix->DefineElement(0, element, 0.01983871336);
464 element = table->GetElement(6);
465 mix->DefineElement(1, element, 0.520088819984);
466 element = table->GetElement(7);
467 mix->DefineElement(2, element, 0.0551367996);
468 element = table->GetElement(8);
469 mix->DefineElement(3, element, 0.157399667056);
470 element = table->GetElement(13);
471 mix->DefineElement(4, element, 0.247536);
473 case 13: //SDD LV+signal cable
474 mix = new TGeoMixture("ITS_LV-CABLE", 5, 2.1035);
475 element = table->GetElement(1);
476 mix->DefineElement(0, element, 0.0082859922);
477 element = table->GetElement(6);
478 mix->DefineElement(1, element, 0.21722436468);
479 element = table->GetElement(7);
480 mix->DefineElement(2, element, 0.023028867);
481 element = table->GetElement(8);
482 mix->DefineElement(3, element, 0.06574077612);
483 element = table->GetElement(13);
484 mix->DefineElement(4, element, 0.68572);
486 case 14: //SDD hybrid microcab
487 mix = new TGeoMixture("ITS_HYB-CAB", 5, 2.0502);
488 element = table->GetElement(1);
489 mix->DefineElement(0, element, 0.00926228815);
490 element = table->GetElement(6);
491 mix->DefineElement(1, element, 0.24281879711);
492 element = table->GetElement(7);
493 mix->DefineElement(2, element, 0.02574224025);
494 element = table->GetElement(8);
495 mix->DefineElement(3, element, 0.07348667449);
496 element = table->GetElement(13);
497 mix->DefineElement(4, element, 0.64869);
499 case 15: //SDD anode microcab
500 mix = new TGeoMixture("ITS_ANOD-CAB", 5, 1.7854);
501 element = table->GetElement(1);
502 mix->DefineElement(0, element, 0.0128595919215);
503 element = table->GetElement(6);
504 mix->DefineElement(1, element, 0.392653705471);
505 element = table->GetElement(7);
506 mix->DefineElement(2, element, 0.041626868025);
507 element = table->GetElement(8);
508 mix->DefineElement(3, element, 0.118832707289);
509 element = table->GetElement(13);
510 mix->DefineElement(4, element, 0.431909);
512 case 16: // inox/alum
513 mix = new TGeoMixture("ITS_INOX-AL", 5, 3.0705);
514 element = table->GetElement(13);
515 mix->DefineElement(0, element, 0.816164);
516 element = table->GetElement(14);
517 mix->DefineElement(1, element, 0.000919182);
518 element = table->GetElement(24);
519 mix->DefineElement(2, element, 0.0330906);
520 element = table->GetElement(26);
521 mix->DefineElement(3, element, 0.131443);
522 element = table->GetElement(28);
523 mix->DefineElement(4, element, 0.0183836);
525 mix = new TGeoMixture("TPC_MYLAR", 3, 1.39);
526 element = table->GetElement(1);
527 mix->DefineElement(0, element, 0.0416667);
528 element = table->GetElement(6);
529 mix->DefineElement(1, element, 0.625);
530 element = table->GetElement(8);
531 mix->DefineElement(2, element, 0.333333);
533 case 18: // SPDBUS(AL+KPT+EPOX) - unknown composition
534 mix = new TGeoMixture("ITS_SPDBUS", 1, 1.906);
535 element = table->GetElement(9);
536 mix->DefineElement(0, element, 1.);
539 case 19: // SDD/SSD rings - unknown composition
540 mix = new TGeoMixture("ITS_SDDRINGS", 1, 1.8097);
541 element = table->GetElement(6);
542 mix->DefineElement(0, element, 1.);
545 case 20: // SPD end ladder - unknown composition
546 mix = new TGeoMixture("ITS_SPDEL", 1, 3.6374);
547 element = table->GetElement(22);
548 mix->DefineElement(0, element, 1.);
551 case 21: // SDD end ladder - unknown composition
552 mix = new TGeoMixture("ITS_SDDEL", 1, 0.3824);
553 element = table->GetElement(30);
554 mix->DefineElement(0, element, 1.);
557 case 22: // SSD end ladder - unknown composition
558 mix = new TGeoMixture("ITS_SSDEL", 1, 0.68);
559 element = table->GetElement(16);
560 mix->DefineElement(0, element, 1.);
565 printf("Patched with mixture %s\n", mix->GetName());
569 //_____________________________________________________________________________
570 void TFlukaMCGeometry::CreateFlukaMatFile(const char *fname)
572 // ==== from FLUGG ====
573 // NAMES OF ELEMENTS AND COMPOUNDS: the names must be written in upper case,
574 // according to the fluka standard. In addition,. they must be equal to the
575 // names of the fluka materials - see fluka manual - in order that the
576 // program load the right cross sections, and equal to the names included in
577 // the .pemf. Otherwise the user must define the LOW-MAT CARDS, and make his
578 // own .pemf, in order to get the right cross sections loaded in memory.
582 gGeoManager->Export("flgeom.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 // element = table->GetElement(10);
603 // element->SetTitle("ARGON"); // NEON not in neutron xsec table
604 Int_t nelements = table->GetNelements();
605 TList *matlist = gGeoManager->GetListOfMaterials();
606 // TList *medlist = gGeoManager->GetListOfMedia();
607 // Int_t nmed = medlist->GetSize();
609 Int_t nmater = matlist->GetSize();
612 TGeoMixture *mix = 0;
615 // Create all needed elements
616 for (Int_t i=1; i<nelements; i++) {
617 element = table->GetElement(i);
618 // skip elements which are not defined
619 if (!element->IsUsed() && !element->IsDefined()) continue;
620 matname = element->GetTitle();
621 ToFlukaString(matname);
624 mat = new TGeoMaterial(matname, element->A(), element->Z(), rho);
625 mat->SetIndex(nfmater+3);
628 objstr = new TObjString(matname.Data());
629 fMatNames->Add(objstr);
635 // Adjust material names and add them to FLUKA list
636 for (i=0; i<nmater; i++) {
637 mat = (TGeoMaterial*)matlist->At(i);
638 if (!mat->IsUsed()) continue;
641 rho = mat->GetDensity();
642 if (mat->GetZ()<0.001) {
643 mat->SetIndex(2); // vacuum, built-in inside FLUKA
646 matname = mat->GetName();
647 FlukaMatName(matname);
649 mat->SetIndex(nfmater+3);
650 objstr = new TObjString(matname.Data());
652 fMatNames->Add(objstr);
656 // Dump all elements with MATERIAL cards
657 for (i=0; i<nfmater; i++) {
658 mat = (TGeoMaterial*)fMatList->At(i);
659 // mat->SetUsed(kFALSE);
661 out << setw(10) << "MATERIAL ";
662 out.setf(static_cast<std::ios::fmtflags>(0),std::ios::floatfield);
663 objstr = (TObjString*)fMatNames->At(i);
664 matname = objstr->GetString();
667 rho = mat->GetDensity();
668 if (mat->IsMixture()) {
669 out << setw(10) << " ";
670 out << setw(10) << " ";
671 mix = (TGeoMixture*)mat;
673 out << setw(10) << setiosflags(ios::fixed) << setprecision(1) << z;
674 out << setw(10) << setprecision(3) << a;
676 out.setf(static_cast<std::ios::fmtflags>(0),std::ios::floatfield);
677 out << setw(10) << setiosflags(ios::scientific) << setprecision(3) << rho;
678 out.setf(static_cast<std::ios::fmtflags>(0),std::ios::floatfield);
679 out << setw(10) << setiosflags(ios::fixed) << setprecision(1) << Double_t(mat->GetIndex());
680 out << setw(10) << " ";
681 out << setw(10) << " ";
682 out << setw(8) << matname.Data() << endl;
684 // add LOW-MAT card for NEON to associate with ARGON neutron xsec
686 out << setw(10) << "LOW-MAT ";
687 out.setf(static_cast<std::ios::fmtflags>(0),std::ios::floatfield);
688 out << setw(10) << setiosflags(ios::fixed) << setprecision(1) << Double_t(mat->GetIndex());
689 out << setw(10) << setiosflags(ios::fixed) << setprecision(1) << 18.;
690 out << setw(10) << setiosflags(ios::fixed) << setprecision(1) << -2.;
691 out << setw(10) << setiosflags(ios::fixed) << setprecision(1) << 293.;
692 out << setw(10) << " ";
693 out << setw(10) << " ";
694 // out << setw(8) << matname.Data() << endl;
695 out << setw(8) << " " << endl;
698 element = table->GetElement((int)z);
699 TString elename = element->GetTitle();
700 ToFlukaString(elename);
701 if( matname.CompareTo( elename ) != 0 ) {
702 out << setw(10) << "LOW-MAT ";
703 out.setf(static_cast<std::ios::fmtflags>(0),std::ios::floatfield);
704 out << setw(10) << setiosflags(ios::fixed) << setprecision(1) << Double_t(mat->GetIndex());
705 out << setw(10) << setiosflags(ios::fixed) << setprecision(1) << z;
706 out << setw(10) << setiosflags(ios::fixed) << setprecision(1) << " ";
707 out << setw(10) << setiosflags(ios::fixed) << setprecision(1) << " ";
708 out << setw(10) << " ";
709 out << setw(10) << " ";
710 // missing material at Low Energy Cross Section Table
711 if( (int)z==10 || (int)z==21 || (int)z==34 || (int)z==37 || (int)z==39 || (int)z==44 ||
712 (int)z==45 || (int)z==46 || (int)z==52 || (int)z==57 || (int)z==59 || (int)z==60 ||
713 (int)z==61 || (int)z==65 || (int)z==66 || (int)z==67 || (int)z==68 || (int)z==69 ||
714 (int)z==70 || (int)z==71 || (int)z==72 || (int)z==76 || (int)z==77 || (int)z==78 ||
715 (int)z==81 || (int)z==84 || (int)z==85 || (int)z==86 || (int)z==87 || (int)z==88 ||
716 (int)z==89 || (int)z==91 )
717 out << setw(8) << "UNKNOWN " << endl;
719 out << setw(8) << elename.Data() << endl;
720 // out << setw(8) << " " << endl;
726 out << setw(10) << "COMPOUND ";
727 nelem = mix->GetNelements();
728 objstr = (TObjString*)fMatNames->At(i);
729 matname = objstr->GetString();
730 for (j=0; j<nelem; j++) {
731 w = (mix->GetWmixt())[j];
732 if (w<0.00001) w=0.00001;
733 z = (mix->GetZmixt())[j];
734 a = (mix->GetAmixt())[j];
735 idmat = GetElementIndex(Int_t(z));
736 if (!idmat) Error("CreateFlukaMatFile", "element with Z=%f not found", z);
737 out.setf(static_cast<std::ios::fmtflags>(0),std::ios::floatfield);
738 out << setw(10) << setiosflags(ios::fixed) << setprecision(6) << -w;
739 out.setf(static_cast<std::ios::fmtflags>(0),std::ios::floatfield);
740 out << setw(10) << setiosflags(ios::fixed) << setprecision(1) << Double_t(idmat);
742 if (counttothree == 3) {
743 out << matname.Data();
745 if ( (j+1) != nelem) out << setw(10) << "COMPOUND ";
750 for (j=0; j<(3-(nelem%3)); j++)
751 out << setw(10) << " " << setw(10) << " ";
752 out << matname.Data();
756 Int_t nvols = gGeoManager->GetListOfUVolumes()->GetEntriesFast()-1;
758 // Now print the material assignments
759 Double_t flagfield = 0.;
760 printf("#############################################################\n");
761 if (gFluka->IsFieldEnabled()) {
763 printf("Magnetic field enabled\n");
764 } else printf("Magnetic field disabled\n");
765 printf("#############################################################\n");
767 PrintHeader(out, "TGEO MATERIAL ASSIGNMENTS");
768 for (i=1; i<=nvols; i++) {
770 vol = gGeoManager->GetVolume(i);
771 mat = vol->GetMedium()->GetMaterial();
772 idmat = mat->GetIndex();
773 for (Int_t j=0; j<nfmater; j++) {
774 mat = (TGeoMaterial*)fMatList->At(j);
775 if (mat->GetIndex() == idmat) mat->SetUsed(kTRUE);
778 Float_t hasfield = (vol->GetMedium()->GetParam(1) > 0) ? flagfield : 0.;
779 out << "* Assigning material: " << vol->GetMedium()->GetMaterial()->GetName() << " to Volume: " << vol->GetName();
782 out << setw(10) << "ASSIGNMAT ";
783 out.setf(static_cast<std::ios::fmtflags>(0),std::ios::floatfield);
784 out << setw(10) << setiosflags(ios::fixed) << Double_t(idmat);
785 out << setw(10) << setiosflags(ios::fixed) << Double_t(i);
786 out << setw(10) << "0.0";
787 out << setw(10) << "0.0";
788 out << setw(10) << setiosflags(ios::fixed) << hasfield;
789 out << setw(10) << "0.0";
794 fDummyRegion = nvols+1;
795 out << "* Dummy region: " << endl;
796 out << setw(10) << "ASSIGNMAT ";
797 out.setf(static_cast<std::ios::fmtflags>(0),std::ios::floatfield);
798 out << setw(10) << setiosflags(ios::fixed) << idmat;
799 out << setw(10) << setiosflags(ios::fixed) << fDummyRegion;
800 out << setw(10) << "0.0";
801 out << setw(10) << "0.0";
802 out << setw(10) << "0.0";
803 out << setw(10) << "0.0" << endl;
805 fLastMaterial = nfmater+2;
808 void TFlukaMCGeometry::CreatePemfFile()
811 // Steering routine to write and process peg files producing the pemf input
814 Int_t countMatOK = 0;
815 Int_t countElemError = 0;
816 Int_t countNoStern = 0;
817 Int_t countMixError = 0;
819 Int_t countPemfError = 0;
821 TGeoMaterial* mat = 0x0;
824 for (i = fIndmat; i < fLastMaterial - 2; i++) {
825 printf("Write Peg Files %d\n", i);
827 mat = (TGeoMaterial*)fMatList->At(i);
828 if (!mat->IsUsed()) continue;
830 sprintf(number, "%d", i);
831 sname.Append(number);
834 cout << "******************************************************************************" << endl;
835 cout << "******************************************************************************" << endl;
837 WritePegFile(i, &countNoStern, &countElemError, &countMixError, &countGas);
838 sname.Prepend("$FLUPRO/pemf/rpemf peg/");
839 gSystem->Exec(sname.Data());
841 // check if the pemf file was created
842 TString sname = Form("peg/mat%d.pemf", i);
843 ifstream in( sname.Data() );
848 cout << "ERROR Fail to create the pemf file " << sname << endl;
852 cout << "Materials (pemf created) " << countMatOK << endl;
853 cout << "Not Sternheimer par. found " << countNoStern << endl;
854 cout << "Elements with error definitions (Z not integer) " << countElemError << endl;
855 cout << "Mixtures with error definitions (Z not integer) " << countMixError << endl;
856 cout << "Posible Gas (rho < 0.01) " << countGas << endl;
857 // cout << "Posible Gas (without pressure information) " << countGasError << endl;
858 cout << "Pemf files Error " << countPemfError << endl;
859 cout << endl << endl;
861 sname = "cat peg/*.pemf > peg/FlukaVmc.pemf";
862 gSystem->Exec(sname.Data());
863 sname = "mv peg/FlukaVmc.pemf FlukaVmc.pemf";
864 gSystem->Exec(sname.Data());
867 //_____________________________________________________________________________
868 void TFlukaMCGeometry::WritePegFile(Int_t imat, Int_t *NoStern, Int_t *ElemError,
869 Int_t *MixError, Int_t *countGas) const
871 // Write the .peg file for one material
873 TGeoMaterial *mat = (TGeoMaterial*)fMatList->At(imat);
874 TString name = ((TObjString*)fMatNames->At(imat))->GetString();
877 TGeoElement *elem = mat->GetElement();
879 TString sname = "mat";
880 sprintf(number, "%d", imat);
881 sname.Append(number);
882 sname.Append(".peg");
883 sname.Prepend("peg/");
885 out.open(sname.Data(), ios::out);
886 if (!out.good()) return;
887 Double_t dens = mat->GetDensity();
888 TGeoMixture *mix = 0;
891 if (mat->IsMixture()) {
892 mix = (TGeoMixture*)mat;
893 nel = mix->GetNelements();
897 cout << "( Element ) " << name << " Z=" << mat->GetZ() << " Rho " << mat->GetDensity() << endl;
899 Double_t zel = mat->GetZ();
900 if( (zel-Int_t(zel))>0.001 || zel < 1 ) {
901 cout << " ERROR: A Element with not integer Z=" << zel << endl;
907 out << "ELEM" << endl;
908 out << " &INP IRAYL=1, RHO=" << dens << ", " << endl;
910 // check for the Sternheimer parameters
911 Double_t *issb_parm = GetISSB( mat->GetDensity(), 1, &zel, 0 );
912 if( issb_parm[0] > 0 && issb_parm[1] > 0 ) {
913 cout << "Sternheimer parameters found" << endl;
914 out << ", ISSB=1, IEV=" << issb_parm[0] << ", CBAR=" << issb_parm[1]
915 << ", X0=" << issb_parm[2] << "," << endl;
916 out << "X1=" <<issb_parm[3] <<", AFACT="<<issb_parm[4] <<", SK="
917 << issb_parm[5] << ", DELTA0=" << issb_parm[6];
920 cout << "WARNING: Strange element, Sternheimer parameters not found" << endl;
926 out << " GASP=1." << endl;
929 out << " &END" << endl;
930 out << name.Data() << endl;
931 out << elem->GetName() << endl;
936 cout << "( Mixture ) " << name << " Rho " << dens << " nElem " << nel << endl;
938 Double_t *zt = new Double_t[nel];
939 Double_t *wt = new Double_t[nel];
940 for (int j=0; j<nel; j++) {
941 zt[j] = (mix->GetZmixt())[j];
942 wt[j] = (mix->GetWmixt())[j];
943 if( (zt[j]-Int_t(zt[j])) > 0.001 || zt[j] < 1 ) {
944 cout << "ERROR Mixture " << name << " with an element with not integer Z=" << zt[j] << endl;
947 // just continue since the mixtures are not patch,
948 // but the final release should include the return
952 Double_t *issb_parm = GetISSB( mat->GetDensity(), nel, zt, wt );
953 out << "MIXT" << endl;
954 out << " &INP IRAYL=1, NE=" << nel << ", RHOZ=" << wt[0] << ",";
955 line = Form(" &INP IRAYL=1, NE=%d RHOZ=%g", nel, wt[0]);
956 for(int j=1; j<nel; j++) {
957 out << " " << wt[j] << ",";
958 line += Form(" %g,", wt[j] );
959 if( line.Length() > 60 ) { out << endl; line = ""; }
961 out << " RHO=" << mat->GetDensity() << ", ";
962 line += Form(" RHO=%g, ", mat->GetDensity());
963 if( line.Length() > 60 ) { out << endl; line = ""; }
965 if( issb_parm[0] > 0 && issb_parm[1] > 0 ) {
966 cout << "Sternheimer parameters found" << endl;
967 out << " ISSB=1, IEV=" << issb_parm[0] << ",";
968 line += Form(" ISSB=1, IEV=%g,", issb_parm[0]);
969 if( line.Length() > 60 ) { out << endl; line = ""; }
970 out << " CBAR=" << issb_parm[1] << ",";
971 line += Form(" CBAR=%g,",issb_parm[1]);
972 if( line.Length() > 60 ) { out << endl; line = ""; }
973 out << " X0=" << issb_parm[2] << ",";
974 line += Form(" X0=%g,", issb_parm[2]);
975 if( line.Length() > 60 ) { out << endl; line = ""; }
976 out << " X1=" << issb_parm[3] << ",";
977 line += Form(" X1=%g,", issb_parm[3]);
978 if( line.Length() > 60 ) { out << endl; line = ""; }
979 out << " AFACT="<< issb_parm[4] << ",";
980 line += Form(" AFACT=%g,", issb_parm[4]);
981 if( line.Length() > 60 ) { out << endl; line = ""; }
982 out << " SK=" << issb_parm[5] << ",";
983 line += Form(" SK=%g,", issb_parm[5]);
984 if( line.Length() > 60 ) { out << endl; line = ""; }
987 cout << "Sternheimer parameters not found" << endl;
993 out << " GASP=1." << endl;
996 out << " &END" << endl;
997 out << name.Data() << endl;
998 for (i=0; i<nel; i++) {
999 elem = mix->GetElement(i);
1000 line = elem->GetName();
1001 if (line.Length()==1) line.Append(" ");
1002 out << line.Data() << " ";
1010 Double_t ue = 3000000.; // [MeV]
1011 Double_t up = 3000000.; // [MeV]
1016 TObjArray* cutList = ((TFluka*) gMC)->GetListOfUserConfigs();
1017 TIter next(cutList);
1018 TFlukaConfigOption* proc;
1020 while((proc = (TFlukaConfigOption*)next()))
1022 if (proc->Medium() == mat->GetIndex()) {
1023 ap = proc->Cut(kCUTGAM);
1024 ae = proc->Cut(kCUTELE);
1025 if (ap == -1.) ap = TFlukaConfigOption::DefaultCut(kCUTGAM);
1026 if (ae == -1.) ae = TFlukaConfigOption::DefaultCut(kCUTELE);
1031 if (ap == -1.) ap = TFlukaConfigOption::DefaultCut(kCUTGAM);
1032 if (ae == -1.) ae = TFlukaConfigOption::DefaultCut(kCUTELE);
1034 ap *= 1000.; // [MeV]
1035 ae = (ae + 0.00051099906) * 1000.; // [MeV]
1037 out << "ENER" << endl;
1038 out << " $INP AE=" << ae << ", UE=" << ue <<", AP=" << ap << ", UP=" << up << " $END" << endl;
1039 out << "PWLF" << endl;
1040 out << " $INP NALE=300, NALG=400, NALR=100 $END" << endl;
1041 out << "DECK" << endl;
1042 out << " $INP $END" << endl;
1043 out << "TEST" << endl;
1044 out << " $INP $END" << endl;
1048 Double_t * TFlukaMCGeometry::GetISSB(Double_t rho, Int_t nElem, Double_t *zelem, Double_t *welem ) const
1050 // Read the density effect parameters
1051 // from R.M. Sternheimer et al. Atomic Data
1052 // and Nuclear Data Tables, Vol. 30 No. 2
1054 // return the parameters if the element/mixture match with one of the list
1055 // otherwise returns the parameters set to 0
1057 struct sternheimerData {
1058 TString longname; // element/mixture name
1059 Int_t nelems; // number of constituents N
1060 Int_t Z[20]; //[nelems] Z
1061 Double_t wt[20]; //[nelems] weight fraction
1062 Double_t density; // g/cm3
1063 Double_t iev; // Average Ion potential (eV)
1064 // **** Sternheimer parameters ****
1065 Double_t cbar; // CBAR
1068 Double_t afact; // AFACT
1070 Double_t delta0; // DELTA0
1078 static Double_t parameters[7];
1079 memset( parameters, 0, sizeof(Double_t) );
1081 static sternheimerData sternDataArray[300];
1082 static Bool_t isFileRead = kFALSE;
1084 // Read the data file if is needed
1085 if( isFileRead == kFALSE ) {
1086 TString sSternheimerInp = getenv("ALICE_ROOT");
1087 sSternheimerInp +="/TFluka/input/Sternheimer.data";
1089 ifstream in(sSternheimerInp);
1091 in.getline(line, 100);
1092 in.getline(line, 100);
1093 in.getline(line, 100);
1094 in.getline(line, 100);
1095 in.getline(line, 100);
1096 in.getline(line, 100);
1100 while( !in.eof() ) {
1101 in >> shortname >> num >> sternDataArray[is].nelems
1102 >> sternDataArray[is].longname >> formula >> state;
1103 if( in.eof() ) break;
1104 for(int i=0; i<sternDataArray[is].nelems; i++) {
1105 in >> sternDataArray[is].Z[i] >> sternDataArray[is].wt[i];
1107 in >> sternDataArray[is].density;
1108 in >> sternDataArray[is].iev;
1109 in >> sternDataArray[is].cbar;
1110 in >> sternDataArray[is].x0;
1111 in >> sternDataArray[is].x1;
1112 in >> sternDataArray[is].afact;
1113 in >> sternDataArray[is].sk;
1114 if( sternDataArray[is].nelems == 1 ) in >> sternDataArray[is].delta0;
1124 // check for elements
1125 if( sternDataArray[is].nelems == 1 && nElem == 1
1126 && sternDataArray[is].Z[0] == Int_t(*zelem)
1127 && TMath::Abs( (sternDataArray[is].density - rho)/sternDataArray[is].density ) < 0.1 ) {
1128 cout << sternDataArray[is].longname << " #elems:" << sternDataArray[is].nelems << " Rho:"
1129 << sternDataArray[is].density << endl;
1130 cout << sternDataArray[is].iev << " "
1131 << sternDataArray[is].cbar << " "
1132 << sternDataArray[is].x0 << " "
1133 << sternDataArray[is].x1 << " "
1134 << sternDataArray[is].afact << " "
1135 << sternDataArray[is].sk << " "
1136 << sternDataArray[is].delta0 << endl;
1138 parameters[0] = sternDataArray[is].iev;
1139 parameters[1] = sternDataArray[is].cbar;
1140 parameters[2] = sternDataArray[is].x0;
1141 parameters[3] = sternDataArray[is].x1;
1142 parameters[4] = sternDataArray[is].afact;
1143 parameters[5] = sternDataArray[is].sk;
1144 parameters[6] = sternDataArray[is].delta0;
1148 // check for mixture
1150 if( sternDataArray[is].nelems > 1 && sternDataArray[is].nelems == nElem ) {
1151 for(int j=0; j<sternDataArray[is].nelems; j++) {
1152 if( sternDataArray[is].Z[j] == Int_t(zelem[j]) &&
1153 TMath::Abs( (sternDataArray[is].wt[j] - welem[j])/sternDataArray[is].wt[j] ) < 0.1 )
1158 if( sternDataArray[is].nelems > 1 &&
1159 TMath::Abs( (sternDataArray[is].density - rho)/sternDataArray[is].density ) < 0.1
1160 && nmatch == sternDataArray[is].nelems ) {
1161 cout << sternDataArray[is].longname << " #elem:" << sternDataArray[is].nelems << " Rho:"
1162 << sternDataArray[is].density << endl;
1163 cout << sternDataArray[is].iev << " "
1164 << sternDataArray[is].cbar << " "
1165 << sternDataArray[is].x0 << " "
1166 << sternDataArray[is].x1 << " "
1167 << sternDataArray[is].afact << " "
1168 << sternDataArray[is].sk << " "
1169 << sternDataArray[is].delta0 << endl;
1171 parameters[0] = sternDataArray[is].iev;
1172 parameters[1] = sternDataArray[is].cbar;
1173 parameters[2] = sternDataArray[is].x0;
1174 parameters[3] = sternDataArray[is].x1;
1175 parameters[4] = sternDataArray[is].afact;
1176 parameters[5] = sternDataArray[is].sk;
1185 //_____________________________________________________________________________
1186 void TFlukaMCGeometry::PrintHeader(ofstream &out, const char *text) const
1188 // Print a FLUKA header.
1189 out << "*\n" << "*\n" << "*\n";
1190 out << "********************* " << text << " *********************\n"
1192 out << "*...+....1....+....2....+....3....+....4....+....5....+....6....+....7..."
1197 //_____________________________________________________________________________
1198 Int_t TFlukaMCGeometry::RegionId() const
1200 // Returns current region id <-> TGeo node id
1201 if (gGeoManager->IsOutside()) return 0;
1202 return gGeoManager->GetCurrentNode()->GetUniqueID();
1205 //_____________________________________________________________________________
1206 Int_t TFlukaMCGeometry::GetElementIndex(Int_t z) const
1208 // Get index of a material having a given Z element.
1209 TIter next(fMatList);
1212 while ((mat=(TGeoMaterial*)next())) {
1213 if (mat->IsMixture()) continue;
1214 if (mat->GetElement()->Z() == z) return mat->GetIndex();
1219 //_____________________________________________________________________________
1220 void TFlukaMCGeometry::SetMreg(Int_t mreg, Int_t lttc)
1222 // Update if needed next history;
1223 // if (gFluka->GetDummyBoundary()==2) {
1224 // gGeoManager->CdNode(fNextLattice-1);
1227 if (lttc == TFlukaMCGeometry::kLttcOutside) {
1228 fCurrentRegion = NofVolumes()+2;
1229 fCurrentLattice = lttc;
1230 gGeoManager->CdTop();
1231 gGeoManager->SetOutside(kTRUE);
1233 if (lttc == TFlukaMCGeometry::kLttcVirtual) return;
1235 Error("TFlukaMCGeometry::SetMreg","Invalide lattice %i",lttc);
1238 fCurrentRegion = mreg;
1239 fCurrentLattice = lttc;
1241 Int_t crtlttc = gGeoManager->GetCurrentNodeId()+1;
1242 if (crtlttc == lttc) return;
1243 gGeoManager->CdNode(lttc-1);
1246 //_____________________________________________________________________________
1247 void TFlukaMCGeometry::SetCurrentRegion(Int_t mreg, Int_t latt)
1249 // Set index/history for next entered region
1250 fCurrentRegion = mreg;
1251 fCurrentLattice = latt;
1254 //_____________________________________________________________________________
1255 void TFlukaMCGeometry::SetNextRegion(Int_t mreg, Int_t latt)
1257 // Set index/history for next entered region
1259 fNextLattice = latt;
1262 //_____________________________________________________________________________
1263 void TFlukaMCGeometry::ToFlukaString(TString &str) const
1265 // ToFlukaString converts an string to something usefull in FLUKA:
1266 // * Capital letters
1268 // * Replace ' ' by '_'
1269 if (str.Length()<8) {
1274 for (ilast=7; ilast>0; ilast--) if (str(ilast)!=' ') break;
1276 for (Int_t pos=0; pos<ilast; pos++)
1277 if (str(pos)==' ') str.Replace(pos,1,"_",1);
1281 //_____________________________________________________________________________
1282 void TFlukaMCGeometry::FlukaMatName(TString &str) const
1284 // Strip the detector name
1286 TObjArray * tokens = str.Tokenize("_");
1287 Int_t ntok = tokens->GetEntries();
1289 TString head = ((TObjString*) tokens->At(0))->GetString();
1290 Int_t nhead = head.Length();
1291 str = str.Remove(0, nhead + 1);
1296 // Convert a name to upper case 8 chars.
1299 for (ilast=7; ilast>0; ilast--) if (str(ilast)!=' ') break;
1300 if (ilast>5) ilast = 5;
1302 TIter next(fMatNames);
1306 while ((objstr=(TObjString*)next())) {
1307 matname = objstr->GetString();
1308 if (matname == str) {
1312 sprintf(&number[1], "%d", index);
1313 } else if (index<100) {
1314 sprintf(number, "%d", index);
1316 Error("FlukaMatName", "Too many materials %s", str.Data());
1319 str.Replace(ilast+1, 2, number);
1325 //______________________________________________________________________________
1326 void TFlukaMCGeometry::Vname(const char *name, char *vname) const
1329 // convert name to upper case. Make vname at least 4 chars
1331 Int_t l = strlen(name);
1334 for (i=0;i<l;i++) vname[i] = toupper(name[i]);
1335 for (i=l;i<4;i++) vname[i] = ' ';
1341 // FLUKA GEOMETRY WRAPPERS - to replace FLUGG wrappers
1343 //_____________________________________________________________________________
1344 Int_t idnrwr(const Int_t & /*nreg*/, const Int_t & /*mlat*/)
1347 // Wrapper for setting DNEAR option on fluka side. Must return 0
1348 // if user doesn't want Fluka to use DNEAR to compute the
1349 // step (the same effect is obtained with the GLOBAL (WHAT(3)=-1)
1350 // card in fluka input), returns 1 if user wants Fluka always to
1351 // use DNEAR (in this case, be sure that GEANT4 DNEAR is unique,
1352 // coming from all directions!!!)
1353 if (gMCGeom->IsDebugging()) printf("========== Dummy IDNRWR\n");
1357 //_____________________________________________________________________________
1358 void g1wr(Double_t &pSx, Double_t &pSy, Double_t &pSz,
1359 Double_t *pV, Int_t &oldReg , const Int_t &oldLttc, Double_t &propStep,
1360 Int_t &/*nascFlag*/, Double_t &retStep, Int_t &newReg,
1361 Double_t &saf, Int_t &newLttc, Int_t <tcFlag,
1362 Double_t *sLt, Int_t *jrLt)
1365 // Initialize FLUKa point and direction;
1371 Int_t olttc = oldLttc;
1373 gGeoManager->FindNode(pSx,pSy,pSz);
1374 olttc = gGeoManager->GetCurrentNodeId()+1;
1377 gMCGeom->GetCurrentRegion(ccreg,cclat);
1378 Bool_t crossed = (ccreg==oldReg && cclat==oldLttc)?kFALSE:kTRUE;
1379 gMCGeom->SetCurrentRegion(oldReg, olttc);
1380 // Initialize default return values
1382 jrLt[lttcFlag] = olttc;
1383 sLt[lttcFlag] = propStep;
1384 jrLt[lttcFlag+1] = -1;
1385 sLt[lttcFlag+1] = 0.;
1388 Bool_t crossedDummy = (oldReg == gFluka->GetDummyRegion())?kTRUE:kFALSE;
1389 Int_t curLttc, curReg;
1391 // FLUKA crossed the dummy boundary - update new region/history
1392 retStep = TGeoShape::Tolerance();
1394 gMCGeom->GetNextRegion(newReg, newLttc);
1395 gMCGeom->SetMreg(newReg, newLttc);
1396 sLt[lttcFlag] = TGeoShape::Tolerance(); // null step in current region
1398 jrLt[lttcFlag] = newLttc;
1399 sLt[lttcFlag] = TGeoShape::Tolerance(); // null step in next region
1400 jrLt[lttcFlag+1] = -1;
1401 sLt[lttcFlag+1] = 0.; // null step in next region;
1405 // Reset outside flag
1406 gGeoManager->SetOutside(kFALSE);
1408 curLttc = gGeoManager->GetCurrentNodeId()+1;
1409 curReg = gGeoManager->GetCurrentVolume()->GetNumber();
1410 if (olttc != curLttc) {
1411 // FLUKA crossed the boundary : we trust that the given point is really there,
1412 // so we just update TGeo state
1413 gGeoManager->CdNode(olttc-1);
1414 curLttc = gGeoManager->GetCurrentNodeId()+1;
1415 curReg = gGeoManager->GetCurrentVolume()->GetNumber();
1417 // Now the current TGeo state reflects the FLUKA state
1419 gGeoManager->SetCurrentPoint(pSx, pSy, pSz);
1420 gGeoManager->SetCurrentDirection(pV);
1423 gGeoManager->FindNextBoundaryAndStep(propStep);
1426 gGeoManager->FindNextBoundaryAndStep(propStep, kTRUE);
1427 saf = gGeoManager->GetSafeDistance();
1432 Double_t snext = gGeoManager->GetStep();
1434 if (snext<=0.0) snext = TGeoShape::Tolerance();
1437 NORLAT.distn = snext;
1438 NORLAT.xn[0] += snext*pV[0];
1439 NORLAT.xn[1] += snext*pV[1];
1440 NORLAT.xn[2] += snext*pV[2];
1441 if (!gGeoManager->IsOnBoundary()) {
1442 // Next boundary further than proposed step, which is approved
1443 if (saf>propStep) saf = propStep;
1445 sLt[lttcFlag] = propStep;
1448 if (saf>snext) saf = snext; // Safety should be less than the proposed step if a boundary will be crossed
1449 gGeoManager->SetCurrentPoint(pSx,pSy,pSz);
1450 newLttc = (gGeoManager->IsOutside())?(TFlukaMCGeometry::kLttcOutside):gGeoManager->GetCurrentNodeId()+1;
1451 newReg = (gGeoManager->IsOutside())?(gMCGeom->NofVolumes()+2):gGeoManager->GetCurrentVolume()->GetNumber();
1452 if (gMCGeom->IsDebugging()) printf(" newReg=%i newLttc=%i\n", newReg, newLttc);
1454 // We really crossed the boundary, but is it the same region ?
1455 gMCGeom->SetNextRegion(newReg, newLttc);
1457 Int_t pid = TRACKR.jtrack;
1458 if (newReg==oldReg && newLttc!=olttc && pid!=-1) {
1459 // Virtual boundary between replicants
1460 newReg = gFluka->GetDummyRegion();
1461 newLttc = TFlukaMCGeometry::kLttcVirtual;
1465 sLt[lttcFlag] = snext;
1467 jrLt[lttcFlag] = newLttc;
1468 sLt[lttcFlag] = snext;
1469 jrLt[lttcFlag+1] = -1;
1470 sLt[lttcFlag+1] = 0.;
1471 gGeoManager->SetOutside(kFALSE);
1472 gGeoManager->CdNode(olttc-1);
1473 if (gMCGeom->IsDebugging()) {
1474 printf("=> snext=%g safe=%g\n", snext, saf);
1475 for (Int_t i=0; i<lttcFlag+1; i++) printf(" jrLt[%i]=%i sLt[%i]=%g\n", i,jrLt[i],i,sLt[i]);
1479 //_____________________________________________________________________________
1482 if (gMCGeom->IsDebugging()) printf("========== Dummy G1RTWR\n");
1485 //_____________________________________________________________________________
1486 void conhwr(Int_t & /*intHist*/, Int_t * /*incrCount*/)
1488 if (gMCGeom->IsDebugging()) printf("========== Dummy CONHWR\n");
1491 //_____________________________________________________________________________
1492 void inihwr(Int_t &intHist)
1494 if (gMCGeom->IsDebugging()) printf("========== Inside INIHWR -> reinitializing history: %i\n", intHist);
1495 if (gGeoManager->IsOutside()) gGeoManager->CdTop();
1497 // printf("=== wrong history number\n");
1500 if (intHist==0) gGeoManager->CdTop();
1501 else gGeoManager->CdNode(intHist-1);
1502 if (gMCGeom->IsDebugging()) {
1503 printf(" --- current path: %s\n", gGeoManager->GetPath());
1504 printf("<= INIHWR\n");
1508 //_____________________________________________________________________________
1509 void jomiwr(const Int_t & /*nge*/, const Int_t & /*lin*/, const Int_t & /*lou*/,
1512 // Geometry initialization wrapper called by FLUKAM. Provides to FLUKA the
1513 // number of regions (volumes in TGeo)
1514 // build application geometry
1515 if (gMCGeom->IsDebugging()) printf("========== Inside JOMIWR\n");
1516 flukaReg = gGeoManager->GetListOfUVolumes()->GetEntriesFast()+1;
1517 if (gMCGeom->IsDebugging()) printf("<= JOMIWR: last region=%i\n", flukaReg);
1520 //_____________________________________________________________________________
1521 void lkdbwr(Double_t &pSx, Double_t &pSy, Double_t &pSz,
1522 Double_t *pV, const Int_t &oldReg, const Int_t &oldLttc,
1523 Int_t &flagErr, Int_t &newReg, Int_t &newLttc)
1525 if (gMCGeom->IsDebugging()) {
1526 printf("========== Inside LKDBWR (%f, %f, %f)\n",pSx, pSy, pSz);
1527 printf(" in: pV=(%f, %f, %f)\n", pV[0], pV[1], pV[2]);
1528 printf(" in: oldReg=%i oldLttc=%i\n", oldReg, oldLttc);
1530 lkwr(pSx,pSy,pSz,pV,oldReg,oldLttc,flagErr,newReg,newLttc);
1533 //_____________________________________________________________________________
1534 void lkfxwr(Double_t &pSx, Double_t &pSy, Double_t &pSz,
1535 Double_t *pV, const Int_t &oldReg, const Int_t &oldLttc,
1536 Int_t &flagErr, Int_t &newReg, Int_t &newLttc)
1538 if (gMCGeom->IsDebugging()) {
1539 printf("========== Inside LKFXWR (%f, %f, %f)\n",pSx, pSy, pSz);
1540 printf(" in: pV=(%f, %f, %f)\n", pV[0], pV[1], pV[2]);
1541 printf(" in: oldReg=%i oldLttc=%i\n", oldReg, oldLttc);
1543 lkwr(pSx,pSy,pSz,pV,oldReg,oldLttc,flagErr,newReg,newLttc);
1546 //_____________________________________________________________________________
1547 void lkmgwr(Double_t &pSx, Double_t &pSy, Double_t &pSz,
1548 Double_t *pV, const Int_t &oldReg, const Int_t &oldLttc,
1549 Int_t &flagErr, Int_t &newReg, Int_t &newLttc)
1551 if (gMCGeom->IsDebugging()) {
1552 printf("========== Inside LKMGWR (%f, %f, %f)\n",pSx, pSy, pSz);
1553 printf(" in: pV=(%f, %f, %f)\n", pV[0], pV[1], pV[2]);
1554 printf(" in: oldReg=%i oldLttc=%i\n", oldReg, oldLttc);
1556 lkwr(pSx,pSy,pSz,pV,oldReg,oldLttc,flagErr,newReg,newLttc);
1559 //_____________________________________________________________________________
1560 void lkwr(Double_t &pSx, Double_t &pSy, Double_t &pSz,
1561 Double_t *pV, const Int_t &oldReg, const Int_t &oldLttc,
1562 Int_t &flagErr, Int_t &newReg, Int_t &newLttc)
1564 if (gMCGeom->IsDebugging()) {
1565 printf("========== Inside LKWR (%f, %f, %f)\n",pSx, pSy, pSz);
1566 printf(" in: pV=(%f, %f, %f)\n", pV[0], pV[1], pV[2]);
1567 printf(" in: oldReg=%i oldLttc=%i\n", oldReg, oldLttc);
1570 TGeoNode *node = gGeoManager->FindNode(pSx, pSy, pSz);
1571 if (gGeoManager->IsOutside()) {
1572 newReg = gMCGeom->NofVolumes()+2;
1573 newLttc = TFlukaMCGeometry::kLttcOutside;
1574 gGeoManager->SetOutside(kFALSE);
1575 if (oldLttc>0 && oldLttc<newLttc) gGeoManager->CdNode(oldLttc-1);
1578 gGeoManager->SetOutside(kFALSE);
1579 newReg = node->GetVolume()->GetNumber();
1580 newLttc = gGeoManager->GetCurrentNodeId()+1;
1581 if (oldLttc==TFlukaMCGeometry::kLttcOutside || oldLttc==0) return;
1583 Int_t dummy = gFluka->GetDummyRegion();
1584 if (oldReg==dummy) {
1585 Int_t newreg1, newlttc1;
1586 gMCGeom->GetNextRegion(newreg1, newlttc1);
1587 if (newreg1==newReg && newlttc1==newLttc) {
1589 newLttc = TFlukaMCGeometry::kLttcVirtual;
1594 if (oldReg==newReg && oldLttc!=newLttc) {
1595 newReg = gFluka->GetDummyRegion();
1596 newLttc = TFlukaMCGeometry::kLttcVirtual;
1599 if (gMCGeom->IsDebugging()) {
1600 printf(" LKWR: newReg=%i newLttc=%i\n", newReg, newLttc);
1604 //_____________________________________________________________________________
1605 void nrmlwr(Double_t &pSx, Double_t &pSy, Double_t &pSz,
1606 Double_t &pVx, Double_t &pVy, Double_t &pVz,
1607 Double_t *norml, const Int_t &oldReg,
1608 const Int_t &newReg, Int_t &flagErr)
1610 if (gMCGeom->IsDebugging()) {
1611 printf("========== Inside NRMLWR (%g, %g, %g, %g, %g, %g)\n", pSx,pSy,pSz,pVx,pVy,pVz);
1612 printf(" oldReg=%i, newReg=%i\n", oldReg,newReg);
1614 gGeoManager->SetCurrentPoint(NORLAT.xn[0], NORLAT.xn[1], NORLAT.xn[2]);
1615 gGeoManager->SetCurrentDirection(pVx, pVy, pVz);
1616 Double_t *dnorm = gGeoManager->FindNormalFast();
1619 printf(" ERROR: Cannot compute fast normal\n");
1625 norml[0] = -dnorm[0];
1626 norml[1] = -dnorm[1];
1627 norml[2] = -dnorm[2];
1630 if (gMCGeom->IsDebugging()) {
1631 printf(" normal to boundary: (%g, %g, %g)\n", norml[0], norml[1], norml[2]);
1632 printf("<= NRMLWR\n");
1636 //_____________________________________________________________________________
1637 void rgrpwr(const Int_t & /*flukaReg*/, const Int_t & /*ptrLttc*/, Int_t & /*g4Reg*/,
1638 Int_t * /*indMother*/, Int_t * /*repMother*/, Int_t & /*depthFluka*/)
1640 if (gMCGeom->IsDebugging()) printf("=> Dummy RGRPWR\n");
1643 //_____________________________________________________________________________
1644 Int_t isvhwr(const Int_t &check, const Int_t & intHist)
1647 // Wrapper for saving current navigation history (fCheck=default)
1648 // and returning its pointer. If fCheck=-1 copy of history pointed
1649 // by intHist is made in NavHistWithCount object, and its pointer
1650 // is returned. fCheck=1 and fCheck=2 cases are only in debugging
1651 // version: an array is created by means of FGeometryInit functions
1652 // (but could be a static int * ptrArray = new int[10000] with
1653 // file scope as well) that stores a flag for deleted/undeleted
1654 // histories and at the end of event is checked to verify that
1655 // all saved history objects have been deleted.
1657 // For TGeo, just return the current node ID. No copy need to be made.
1659 if (gMCGeom->IsDebugging()) printf("=> Inside ISVHWR\n");
1660 if (check<0) return intHist;
1661 Int_t histInt = gGeoManager->GetCurrentNodeId()+1;
1662 if (gMCGeom->IsDebugging()) printf("<= ISVHWR: history is: %i in: %s\n", histInt, gGeoManager->GetPath());