Check first if a material (element) has been predefined by FLUKA.
[u/mrichter/AliRoot.git] / TFluka / TFlukaMCGeometry.cxx
CommitLineData
e2e989f9 1/**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3 * *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
6 * *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
15
16// $Id$
c1a0a2f0 17
18// Class TFlukaMCGeometry
19// --------------------
20// Implementation of the TVirtualMCGeometry interface
21// for defining and using TGeo geometry with FLUKA.
22// This allows the FLUKA MonteCarlo to run with the TGeo
23// geometrical modeller
8495a208 24// Author: Andrei Gheata 10/07/2003
25
c1a0a2f0 26#include "Riostream.h"
b16b69d3 27#include "TList.h"
c1a0a2f0 28#include "TCallf77.h"
89125dc7 29#include "TFluka.h"
e44ab76e 30#include "TSystem.h"
8495a208 31#include "TFlukaMCGeometry.h"
6d184c54 32#include "TFlukaConfigOption.h"
8495a208 33#include "TGeoManager.h"
34#include "TGeoVolume.h"
c1a0a2f0 35#include "TObjString.h"
81f1d030 36#include "Fsourcm.h"
aae6bcdd 37#include "Ftrackr.h"
4aba9d66 38#include "Fstepsz.h" //(STEPSZ) fluka common
3860051a 39
8495a208 40#ifndef WIN32
41# define idnrwr idnrwr_
42# define g1wr g1wr_
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_
50# define lkwr lkwr_
51# define magfld magfld_
52# define nrmlwr nrmlwr_
53# define rgrpwr rgrpwr_
54# define isvhwr isvhwr_
55
56#else
57
58# define idnrwr IDNRWR
59# define g1wr G1WR
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
67# define lkwr LKWR
68# define magfld MAGFLD
69# define nrmlwr NRMLWR
70# define rgrpwr RGRPWR
71# define isvhwr ISVHWR
72
73#endif
74
75//____________________________________________________________________________
76extern "C"
77{
78 //
79 // Prototypes for FLUKA navigation methods
80 //
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*/,
4aba9d66 85 Double_t & /*saf*/, Int_t & /*newLttc*/, Int_t & /*LttcFlag*/,
8495a208 86 Double_t *s /*Lt*/, Int_t * /*jrLt*/);
87
88 void type_of_call g1rtwr();
4aba9d66 89 void type_of_call conhwr(Int_t & /*intHist*/, Int_t & /*incrCount*/);
8495a208 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*/,
aae6bcdd 95 Int_t & /*flagErr*/, Int_t & /*newReg*/, Int_t & /*newLttc*/);
8495a208 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*/,
aae6bcdd 98 Int_t & /*flagErr*/, Int_t & /*newReg*/, Int_t & /*newLttc*/);
8495a208 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*/,
4aba9d66 101 Int_t & /*flagErr*/, Int_t & /*newReg*/, Int_t & /*newLttc*/);
8495a208 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*/,
4aba9d66 104 Int_t & /*flagErr*/, Int_t & /*newReg*/, Int_t & /*newLttc*/);
efde9b4d 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*/,
4aba9d66 107// Double_t & /*Bmag*/, Int_t & /*reg*/, Int_t & /*idiscflag*/);
8495a208 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*/,
4aba9d66 110 Double_t * /*norml*/, const Int_t & /*oldReg*/,
111 const Int_t & /*newReg*/, Int_t & /*flagErr*/);
8495a208 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*/);
8be76537 115}
8495a208 116
117// TFluka global pointer
89125dc7 118TFluka *gFluka = 0;
119TFlukaMCGeometry *gMCGeom = 0;
6f5667d1 120Int_t gNstep = 0;
8495a208 121
122ClassImp(TFlukaMCGeometry)
123
81f1d030 124TFlukaMCGeometry* TFlukaMCGeometry::fgInstance= NULL;
8495a208 125
126//_____________________________________________________________________________
127TFlukaMCGeometry::TFlukaMCGeometry(const char *name, const char *title)
4aba9d66 128 :TNamed(name, title),
129 fDebug(kFALSE),
130 fLastMaterial(0),
131 fDummyRegion(0),
132 fCurrentRegion(0),
133 fCurrentLattice(0),
134 fNextRegion(0),
135 fNextLattice(0),
136 fRegionList(0),
137 fIndmat(0),
138 fMatList(new TObjArray(256)),
139 fMatNames(new TObjArray(256))
8495a208 140{
141 //
142 // Standard constructor
143 //
89125dc7 144 gFluka = (TFluka*)gMC;
145 gMCGeom = this;
6f5667d1 146 gNstep = 0;
8495a208 147}
148
149//_____________________________________________________________________________
150TFlukaMCGeometry::TFlukaMCGeometry()
4aba9d66 151 :TNamed(),
152 fDebug(kFALSE),
153 fLastMaterial(0),
154 fDummyRegion(0),
155 fCurrentRegion(0),
156 fCurrentLattice(0),
157 fNextRegion(0),
158 fNextLattice(0),
159 fRegionList(0),
160 fIndmat(0),
161 fMatList(0),
162 fMatNames(0)
163
164{
8495a208 165 //
166 // Default constructor
167 //
89125dc7 168 gFluka = (TFluka*)gMC;
169 gMCGeom = this;
6f5667d1 170 gNstep = 0;
8495a208 171}
172
173//_____________________________________________________________________________
174TFlukaMCGeometry::~TFlukaMCGeometry()
175{
176 //
177 // Destructor
178 //
179 fgInstance=0;
b1536e91 180 if (fRegionList) delete [] fRegionList;
89125dc7 181 if (fMatList) delete fMatList;
182 if (fMatNames) {fMatNames->Delete(); delete fMatNames;}
8495a208 183 if (gGeoManager) delete gGeoManager;
184}
185
186//
187// private methods
188//
8495a208 189
190//_____________________________________________________________________________
191Double_t* TFlukaMCGeometry::CreateDoubleArray(Float_t* array, Int_t size) const
192{
193// Converts Float_t* array to Double_t*,
194// !! The new array has to be deleted by user.
195// ---
196
197 Double_t* doubleArray;
198 if (size>0) {
199 doubleArray = new Double_t[size];
200 for (Int_t i=0; i<size; i++) doubleArray[i] = array[i];
201 }
202 else {
203 //doubleArray = 0;
204 doubleArray = new Double_t[1];
205 }
206 return doubleArray;
207}
208//
209// public methods
8495a208 210
8495a208 211
212//_____________________________________________________________________________
8495a208 213Int_t TFlukaMCGeometry::GetMedium() const
214{
215// Get current medium number
216 Int_t imed = 0;
217 TGeoNode *node = gGeoManager->GetCurrentNode();
218 if (!node) imed = gGeoManager->GetTopNode()->GetVolume()->GetMedium()->GetId();
219 else imed = node->GetVolume()->GetMedium()->GetId();
2bc4c610 220 if (fDebug) printf("GetMedium=%i\n", imed);
8495a208 221 return imed;
222}
223
224//_____________________________________________________________________________
0bc73b6c 225Int_t TFlukaMCGeometry::GetFlukaMaterial(Int_t imed) const
226{
227// Returns FLUKA material index for medium IMED
81f1d030 228 TGeoMedium *med = (TGeoMedium*)gGeoManager->GetListOfMedia()->At(imed-1);
0bc73b6c 229 if (!med) {
230 Error("GetFlukaMaterial", "MEDIUM %i nor found", imed);
231 return -1;
232 }
81f1d030 233 TGeoMaterial* mat = med->GetMaterial();
234 if (!mat->IsUsed()) return -1;
0bc73b6c 235 Int_t imatfl = med->GetMaterial()->GetIndex();
236 return imatfl;
237}
238
239//_____________________________________________________________________________
b1536e91 240Int_t *TFlukaMCGeometry::GetRegionList(Int_t imed, Int_t &nreg)
241{
242// Get an ordered list of regions matching a given medium number
243 nreg = 0;
244 if (!fRegionList) fRegionList = new Int_t[NofVolumes()+1];
245 TIter next(gGeoManager->GetListOfUVolumes());
246 TGeoVolume *vol;
247 Int_t imedium, ireg;
248 while ((vol = (TGeoVolume*)next())) {
e5d91a32 249 TGeoMedium* med;
250 if ((med = vol->GetMedium()) == 0) continue;
251 imedium = med->GetId();
252 if (imedium == imed) {
4aba9d66 253 ireg = vol->GetNumber();
254 fRegionList[nreg++] = ireg;
e5d91a32 255 }
b1536e91 256 }
257 return fRegionList;
4aba9d66 258}
b1536e91 259
260//_____________________________________________________________________________
261Int_t *TFlukaMCGeometry::GetMaterialList(Int_t imat, Int_t &nreg)
262{
263// Get an ordered list of regions matching a given medium number
264 nreg = 0;
265 if (!fRegionList) fRegionList = new Int_t[NofVolumes()+1];
266 TIter next(gGeoManager->GetListOfUVolumes());
267 TGeoVolume *vol;
268 Int_t imaterial, ireg;
269 while ((vol = (TGeoVolume*)next())) {
e5d91a32 270 TGeoMedium* med;
271 if ((med = vol->GetMedium()) == 0) continue;
272 imaterial = med->GetMaterial()->GetIndex();
273 if (imaterial == imat) {
4aba9d66 274 ireg = vol->GetNumber();
275 fRegionList[nreg++] = ireg;
e5d91a32 276 }
b1536e91 277 }
278 return fRegionList;
4aba9d66 279}
8495a208 280
281//_____________________________________________________________________________
8495a208 282Int_t TFlukaMCGeometry::NofVolumes() const
283{
284 //
285 // Return total number of volumes in the geometry
286 //
287
288 return gGeoManager->GetListOfUVolumes()->GetEntriesFast()-1;
289}
8495a208 290
291//_____________________________________________________________________________
89125dc7 292TGeoMaterial * TFlukaMCGeometry::GetMakeWrongMaterial(Double_t z)
293{
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};
298
299 Int_t ind;
300 Double_t dz;
301 for (ind=0; ind<23; ind++) {
302 dz = TMath::Abs(z-kz[ind]);
303 if (dz<1E-4) break;
304 }
305 if (ind>22) {
306 printf("Cannot patch material with Z=%g\n", z);
307 return 0;
308 }
309 TGeoMixture *mix = 0;
310 TGeoElement *element;
81f1d030 311 TGeoElementTable *table = gGeoManager->GetElementTable();
89125dc7 312 switch (ind) {
313 case 0: // AIR
81f1d030 314 mix = new TGeoMixture("TPC_AIR", 4, 0.001205);
89125dc7 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);
323 break;
324 case 1: //SDD SI CHIP
81f1d030 325 mix = new TGeoMixture("ITS_SDD_SI", 6, 2.4485);
89125dc7 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);
338 break;
339 case 2: // WATER
81f1d030 340 mix = new TGeoMixture("ITS_WATER", 2, 1.0);
89125dc7 341 element = table->GetElement(1);
342 mix->DefineElement(0, element, 0.111898344);
343 element = table->GetElement(8);
344 mix->DefineElement(1, element, 0.888101656);
345 break;
346 case 3: // CERAMICS
81f1d030 347 mix = new TGeoMixture("ITS_CERAMICS", 5, 3.6);
89125dc7 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);
358 break;
359 case 4: // EPOXY
81f1d030 360 mix = new TGeoMixture("MUON_G10FR4", 4, 1.8);
89125dc7 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);
369 break;
370 case 5: // EPOXY
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);
380 break;
381 case 6: // KAPTON
81f1d030 382 mix = new TGeoMixture("ITS_KAPTON", 4, 1.3);
89125dc7 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);
391 break;
392 case 7: // INOX
81f1d030 393 mix = new TGeoMixture("ITS_INOX", 9, 7.9);
89125dc7 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);
412 break;
413 case 8: // ROHACELL
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);
423 break;
424 case 9: // SDD-C-AL
81f1d030 425 mix = new TGeoMixture("ITS_SDD-C-AL", 5, 1.9837);
89125dc7 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);
436 break;
437 case 10: // X7R-CAP
81f1d030 438 mix = new TGeoMixture("ITS_X7R-CAP", 7, 6.72);
89125dc7 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);
453 break;
454 case 11: // SDD ruby sph. Al2O3
81f1d030 455 mix = new TGeoMixture("ITS_AL2O3", 2, 3.97);
89125dc7 456 element = table->GetElement(8);
457 mix->DefineElement(0, element, 0.5293);
458 element = table->GetElement(13);
459 mix->DefineElement(1, element, 0.4707);
460 break;
461 case 12: // SDD HV microcable
81f1d030 462 mix = new TGeoMixture("ITS_HV-CABLE", 5, 1.6087);
89125dc7 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);
473 break;
474 case 13: //SDD LV+signal cable
81f1d030 475 mix = new TGeoMixture("ITS_LV-CABLE", 5, 2.1035);
89125dc7 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);
486 break;
487 case 14: //SDD hybrid microcab
81f1d030 488 mix = new TGeoMixture("ITS_HYB-CAB", 5, 2.0502);
89125dc7 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);
499 break;
500 case 15: //SDD anode microcab
81f1d030 501 mix = new TGeoMixture("ITS_ANOD-CAB", 5, 1.7854);
89125dc7 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);
512 break;
513 case 16: // inox/alum
81f1d030 514 mix = new TGeoMixture("ITS_INOX-AL", 5, 3.0705);
89125dc7 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);
525 case 17: // MYLAR
81f1d030 526 mix = new TGeoMixture("TPC_MYLAR", 3, 1.39);
89125dc7 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);
533 break;
534 case 18: // SPDBUS(AL+KPT+EPOX) - unknown composition
81f1d030 535 mix = new TGeoMixture("ITS_SPDBUS", 1, 1.906);
89125dc7 536 element = table->GetElement(9);
537 mix->DefineElement(0, element, 1.);
95161cd7 538 z = element->Z();
89125dc7 539 break;
540 case 19: // SDD/SSD rings - unknown composition
81f1d030 541 mix = new TGeoMixture("ITS_SDDRINGS", 1, 1.8097);
89125dc7 542 element = table->GetElement(6);
543 mix->DefineElement(0, element, 1.);
95161cd7 544 z = element->Z();
89125dc7 545 break;
546 case 20: // SPD end ladder - unknown composition
81f1d030 547 mix = new TGeoMixture("ITS_SPDEL", 1, 3.6374);
89125dc7 548 element = table->GetElement(22);
549 mix->DefineElement(0, element, 1.);
95161cd7 550 z = element->Z();
89125dc7 551 break;
552 case 21: // SDD end ladder - unknown composition
81f1d030 553 mix = new TGeoMixture("ITS_SDDEL", 1, 0.3824);
89125dc7 554 element = table->GetElement(30);
555 mix->DefineElement(0, element, 1.);
95161cd7 556 z = element->Z();
89125dc7 557 break;
558 case 22: // SSD end ladder - unknown composition
81f1d030 559 mix = new TGeoMixture("ITS_SSDEL", 1, 0.68);
89125dc7 560 element = table->GetElement(16);
561 mix->DefineElement(0, element, 1.);
95161cd7 562 z = element->Z();
89125dc7 563 break;
564 }
565 mix->SetZ(z);
566 printf("Patched with mixture %s\n", mix->GetName());
567 return mix;
8495a208 568}
569
570//_____________________________________________________________________________
efde9b4d 571void TFlukaMCGeometry::CreateFlukaMatFile(const char *fname)
8495a208 572{
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.
9968e86c 580 // Materials defined by FLUKA
8495a208 581 TString sname;
582 gGeoManager->Export("flgeom.root");
583 if (fname) sname = fname;
584 else sname = "flukaMat.inp";
585 ofstream out;
586 out.open(sname.Data(), ios::out);
587 if (!out.good()) {
588 Fatal("CreateFlukaMatFile", "could not open file %s for writing", sname.Data());
589 return;
590 }
591 PrintHeader(out, "MATERIALS AND COMPOUNDS");
592 PrintHeader(out, "MATERIALS");
89125dc7 593 Int_t i,j,idmat;
594 Int_t counttothree, nelem;
595 Double_t a,z,rho, w;
81f1d030 596 TGeoElementTable *table = gGeoManager->GetElementTable();
89125dc7 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 ...
89125dc7 602 Int_t nelements = table->GetNelements();
8495a208 603 TList *matlist = gGeoManager->GetListOfMaterials();
604 TIter next(matlist);
605 Int_t nmater = matlist->GetSize();
9968e86c 606 Int_t nfmater = 0;
607 Int_t newind = 26; // here non predefined materials start
608
89125dc7 609 TGeoMaterial *mat;
8495a208 610 TGeoMixture *mix = 0;
611 TString matname;
89125dc7 612 TObjString *objstr;
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);
620 rho = 0.999;
3860051a 621
89125dc7 622 mat = new TGeoMaterial(matname, element->A(), element->Z(), rho);
9968e86c 623 // Check if the element has been predefined by FLUKA
624 Int_t pmid = GetPredefinedMaterialId(Int_t(element->Z()));
625 if (pmid > 0) {
626 mat->SetIndex(pmid);
627 } else {
628 mat->SetIndex(newind++);
629 }
630
89125dc7 631 mat->SetUsed(kTRUE);
632 fMatList->Add(mat);
633 objstr = new TObjString(matname.Data());
634 fMatNames->Add(objstr);
635 nfmater++;
636 }
6d184c54 637
638 fIndmat = nfmater;
89125dc7 639 // Adjust material names and add them to FLUKA list
8495a208 640 for (i=0; i<nmater; i++) {
641 mat = (TGeoMaterial*)matlist->At(i);
81f1d030 642 if (!mat->IsUsed()) continue;
89125dc7 643 z = mat->GetZ();
644 a = mat->GetA();
645 rho = mat->GetDensity();
646 if (mat->GetZ()<0.001) {
efde9b4d 647 mat->SetIndex(2); // vacuum, built-in inside FLUKA
648 continue;
89125dc7 649 }
650 matname = mat->GetName();
651 FlukaMatName(matname);
81f1d030 652
9968e86c 653 mat->SetIndex(newind++);
89125dc7 654 objstr = new TObjString(matname.Data());
655 fMatList->Add(mat);
656 fMatNames->Add(objstr);
657 nfmater++;
658 }
659
660 // Dump all elements with MATERIAL cards
8495a208 661 for (i=0; i<nfmater; i++) {
89125dc7 662 mat = (TGeoMaterial*)fMatList->At(i);
663// mat->SetUsed(kFALSE);
664 mix = 0;
8495a208 665 out << setw(10) << "MATERIAL ";
666 out.setf(static_cast<std::ios::fmtflags>(0),std::ios::floatfield);
89125dc7 667 objstr = (TObjString*)fMatNames->At(i);
b0853588 668 matname = objstr->GetString();
89125dc7 669 z = mat->GetZ();
670 a = mat->GetA();
671 rho = mat->GetDensity();
8495a208 672 if (mat->IsMixture()) {
673 out << setw(10) << " ";
674 out << setw(10) << " ";
675 mix = (TGeoMixture*)mat;
676 } else {
89125dc7 677 out << setw(10) << setiosflags(ios::fixed) << setprecision(1) << z;
678 out << setw(10) << setprecision(3) << a;
8495a208 679 }
680 out.setf(static_cast<std::ios::fmtflags>(0),std::ios::floatfield);
89125dc7 681 out << setw(10) << setiosflags(ios::scientific) << setprecision(3) << rho;
8495a208 682 out.setf(static_cast<std::ios::fmtflags>(0),std::ios::floatfield);
89125dc7 683 out << setw(10) << setiosflags(ios::fixed) << setprecision(1) << Double_t(mat->GetIndex());
8495a208 684 out << setw(10) << " ";
efde9b4d 685 out << setw(10) << " ";
686 out << setw(8) << matname.Data() << endl;
89125dc7 687 if (!mix) {
688 // add LOW-MAT card for NEON to associate with ARGON neutron xsec
95161cd7 689 if (z==10) {
89125dc7 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;
95161cd7 700 }
701 else {
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;
722 else
723 out << setw(8) << elename.Data() << endl;
724 // out << setw(8) << " " << endl;
725 }
726 }
89125dc7 727 continue;
d11bf3ca 728 }
8495a208 729 counttothree = 0;
730 out << setw(10) << "COMPOUND ";
731 nelem = mix->GetNelements();
89125dc7 732 objstr = (TObjString*)fMatNames->At(i);
8495a208 733 matname = objstr->GetString();
8495a208 734 for (j=0; j<nelem; j++) {
89125dc7 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);
8495a208 741 out.setf(static_cast<std::ios::fmtflags>(0),std::ios::floatfield);
89125dc7 742 out << setw(10) << setiosflags(ios::fixed) << setprecision(6) << -w;
8495a208 743 out.setf(static_cast<std::ios::fmtflags>(0),std::ios::floatfield);
744 out << setw(10) << setiosflags(ios::fixed) << setprecision(1) << Double_t(idmat);
745 counttothree++;
746 if (counttothree == 3) {
747 out << matname.Data();
748 out << endl;
749 if ( (j+1) != nelem) out << setw(10) << "COMPOUND ";
750 counttothree = 0;
89125dc7 751 }
752 }
8495a208 753 if (nelem%3) {
754 for (j=0; j<(3-(nelem%3)); j++)
755 out << setw(10) << " " << setw(10) << " ";
756 out << matname.Data();
757 out << endl;
89125dc7 758 }
759 }
8495a208 760 Int_t nvols = gGeoManager->GetListOfUVolumes()->GetEntriesFast()-1;
761 TGeoVolume *vol;
8495a208 762 // Now print the material assignments
d11bf3ca 763 Double_t flagfield = 0.;
764 printf("#############################################################\n");
89125dc7 765 if (gFluka->IsFieldEnabled()) {
d11bf3ca 766 flagfield = 1.;
767 printf("Magnetic field enabled\n");
768 } else printf("Magnetic field disabled\n");
769 printf("#############################################################\n");
770
8495a208 771 PrintHeader(out, "TGEO MATERIAL ASSIGNMENTS");
772 for (i=1; i<=nvols; i++) {
e5d91a32 773 TGeoMedium* med;
8495a208 774 vol = gGeoManager->GetVolume(i);
e5d91a32 775 if ((med = vol->GetMedium()) == 0) continue;
776 mat = med->GetMaterial();
8495a208 777 idmat = mat->GetIndex();
89125dc7 778 for (Int_t j=0; j<nfmater; j++) {
779 mat = (TGeoMaterial*)fMatList->At(j);
780 if (mat->GetIndex() == idmat) mat->SetUsed(kTRUE);
781 }
3860051a 782
783 Float_t hasfield = (vol->GetMedium()->GetParam(1) > 0) ? flagfield : 0.;
89aa6d28 784 out << "* Assigning material: " << vol->GetMedium()->GetMaterial()->GetName() << " to Volume: " << vol->GetName();
785 out << endl;
786
8495a208 787 out << setw(10) << "ASSIGNMAT ";
788 out.setf(static_cast<std::ios::fmtflags>(0),std::ios::floatfield);
789 out << setw(10) << setiosflags(ios::fixed) << Double_t(idmat);
790 out << setw(10) << setiosflags(ios::fixed) << Double_t(i);
791 out << setw(10) << "0.0";
b0853588 792 out << setw(10) << "0.0";
3860051a 793 out << setw(10) << setiosflags(ios::fixed) << hasfield;
b0853588 794 out << setw(10) << "0.0";
8495a208 795 out << endl;
796 }
aae6bcdd 797 // dummy region
798 idmat = 2; // vacuum
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;
8495a208 809 out.close();
efde9b4d 810 fLastMaterial = nfmater+2;
6d184c54 811}
812
813void TFlukaMCGeometry::CreatePemfFile()
814{
815 //
816 // Steering routine to write and process peg files producing the pemf input
817 //
818 char number[20];
819 Int_t countMatOK = 0;
820 Int_t countElemError = 0;
821 Int_t countNoStern = 0;
822 Int_t countMixError = 0;
823 Int_t countGas = 0;
824 Int_t countPemfError = 0;
825 Int_t i;
826 TGeoMaterial* mat = 0x0;
827 TString sname;
828
829 for (i = fIndmat; i < fLastMaterial - 2; i++) {
4aba9d66 830 printf("Write Peg Files %d\n", i);
831
832 mat = (TGeoMaterial*)fMatList->At(i);
833 if (!mat->IsUsed()) continue;
834 sname = "mat";
835 sprintf(number, "%d", i);
836 sname.Append(number);
837 cout << endl;
838 cout << endl;
839 cout << "******************************************************************************" << endl;
840 cout << "******************************************************************************" << endl;
841 cout << endl;
842 WritePegFile(i, &countNoStern, &countElemError, &countMixError, &countGas);
843 sname.Prepend("$FLUPRO/pemf/rpemf peg/");
844 gSystem->Exec(sname.Data());
845
846 // check if the pemf file was created
847 TString sname = Form("peg/mat%d.pemf", i);
848 ifstream in( sname.Data() );
849 if ( in ) {
850 countMatOK++;
851 in.close();
852 } else {
853 cout << "ERROR Fail to create the pemf file " << sname << endl;
854 countPemfError++;
855 }
6d184c54 856 }
857 cout << "Materials (pemf created) " << countMatOK << endl;
858 cout << "Not Sternheimer par. found " << countNoStern << endl;
859 cout << "Elements with error definitions (Z not integer) " << countElemError << endl;
860 cout << "Mixtures with error definitions (Z not integer) " << countMixError << endl;
861 cout << "Posible Gas (rho < 0.01) " << countGas << endl;
862 // cout << "Posible Gas (without pressure information) " << countGasError << endl;
95161cd7 863 cout << "Pemf files Error " << countPemfError << endl;
6d184c54 864 cout << endl << endl;
865
866 sname = "cat peg/*.pemf > peg/FlukaVmc.pemf";
867 gSystem->Exec(sname.Data());
868 sname = "mv peg/FlukaVmc.pemf FlukaVmc.pemf";
869 gSystem->Exec(sname.Data());
89125dc7 870}
871
872//_____________________________________________________________________________
95161cd7 873void TFlukaMCGeometry::WritePegFile(Int_t imat, Int_t *NoStern, Int_t *ElemError,
874 Int_t *MixError, Int_t *countGas) const
89125dc7 875{
95161cd7 876 // Write the .peg file for one material
877
89125dc7 878 TGeoMaterial *mat = (TGeoMaterial*)fMatList->At(imat);
879 TString name = ((TObjString*)fMatNames->At(imat))->GetString();
880 TString line;
881 char number[20];
882 TGeoElement *elem = mat->GetElement();
883 name = name.Strip();
884 TString sname = "mat";
885 sprintf(number, "%d", imat);
886 sname.Append(number);
887 sname.Append(".peg");
888 sname.Prepend("peg/");
889 ofstream out;
890 out.open(sname.Data(), ios::out);
891 if (!out.good()) return;
892 Double_t dens = mat->GetDensity();
893 TGeoMixture *mix = 0;
894 Int_t nel = 1;
895 Int_t i;
896 if (mat->IsMixture()) {
897 mix = (TGeoMixture*)mat;
898 nel = mix->GetNelements();
95161cd7 899 }
900
89125dc7 901 if (nel==1) {
95161cd7 902 cout << "( Element ) " << name << " Z=" << mat->GetZ() << " Rho " << mat->GetDensity() << endl;
903
904 Double_t zel = mat->GetZ();
905 if( (zel-Int_t(zel))>0.001 || zel < 1 ) {
906 cout << " ERROR: A Element with not integer Z=" << zel << endl;
907 cout << endl;
908 (*ElemError)++;
909 return;
910 }
911
89125dc7 912 out << "ELEM" << endl;
913 out << " &INP IRAYL=1, RHO=" << dens << ", " << endl;
95161cd7 914
915 // check for the Sternheimer parameters
916 Double_t *issb_parm = GetISSB( mat->GetDensity(), 1, &zel, 0 );
917 if( issb_parm[0] > 0 && issb_parm[1] > 0 ) {
918 cout << "Sternheimer parameters found" << endl;
919 out << ", ISSB=1, IEV=" << issb_parm[0] << ", CBAR=" << issb_parm[1]
920 << ", X0=" << issb_parm[2] << "," << endl;
921 out << "X1=" <<issb_parm[3] <<", AFACT="<<issb_parm[4] <<", SK="
922 << issb_parm[5] << ", DELTA0=" << issb_parm[6];
923 }
924 else {
925 cout << "WARNING: Strange element, Sternheimer parameters not found" << endl;
926 (*NoStern)++;
927 }
928
929 if (dens<0.01) {
930 (*countGas)++;
931 out << " GASP=1." << endl;
932 }
933
89125dc7 934 out << " &END" << endl;
935 out << name.Data() << endl;
936 out << elem->GetName() << endl;
95161cd7 937
938 }
939 else {
940
941 cout << "( Mixture ) " << name << " Rho " << dens << " nElem " << nel << endl;
942
943 Double_t *zt = new Double_t[nel];
944 Double_t *wt = new Double_t[nel];
945 for (int j=0; j<nel; j++) {
946 zt[j] = (mix->GetZmixt())[j];
947 wt[j] = (mix->GetWmixt())[j];
948 if( (zt[j]-Int_t(zt[j])) > 0.001 || zt[j] < 1 ) {
949 cout << "ERROR Mixture " << name << " with an element with not integer Z=" << zt[j] << endl;
950 cout << endl;
951 (*MixError)++;
952 // just continue since the mixtures are not patch,
953 // but the final release should include the return
954 // return;
955 }
956 }
957 Double_t *issb_parm = GetISSB( mat->GetDensity(), nel, zt, wt );
89125dc7 958 out << "MIXT" << endl;
95161cd7 959 out << " &INP IRAYL=1, NE=" << nel << ", RHOZ=" << wt[0] << ",";
960 line = Form(" &INP IRAYL=1, NE=%d RHOZ=%g", nel, wt[0]);
961 for(int j=1; j<nel; j++) {
962 out << " " << wt[j] << ",";
963 line += Form(" %g,", wt[j] );
964 if( line.Length() > 60 ) { out << endl; line = ""; }
89125dc7 965 }
95161cd7 966 out << " RHO=" << mat->GetDensity() << ", ";
967 line += Form(" RHO=%g, ", mat->GetDensity());
968 if( line.Length() > 60 ) { out << endl; line = ""; }
969
970 if( issb_parm[0] > 0 && issb_parm[1] > 0 ) {
971 cout << "Sternheimer parameters found" << endl;
972 out << " ISSB=1, IEV=" << issb_parm[0] << ",";
973 line += Form(" ISSB=1, IEV=%g,", issb_parm[0]);
974 if( line.Length() > 60 ) { out << endl; line = ""; }
975 out << " CBAR=" << issb_parm[1] << ",";
976 line += Form(" CBAR=%g,",issb_parm[1]);
977 if( line.Length() > 60 ) { out << endl; line = ""; }
978 out << " X0=" << issb_parm[2] << ",";
979 line += Form(" X0=%g,", issb_parm[2]);
980 if( line.Length() > 60 ) { out << endl; line = ""; }
981 out << " X1=" << issb_parm[3] << ",";
982 line += Form(" X1=%g,", issb_parm[3]);
983 if( line.Length() > 60 ) { out << endl; line = ""; }
984 out << " AFACT="<< issb_parm[4] << ",";
985 line += Form(" AFACT=%g,", issb_parm[4]);
986 if( line.Length() > 60 ) { out << endl; line = ""; }
987 out << " SK=" << issb_parm[5] << ",";
988 line += Form(" SK=%g,", issb_parm[5]);
989 if( line.Length() > 60 ) { out << endl; line = ""; }
990 }
991 else {
992 cout << "Sternheimer parameters not found" << endl;
993 (*NoStern)++;
994 }
995
996 if (dens<0.01){
997 (*countGas)++;
998 out << " GASP=1." << endl;
999 }
1000
1001 out << " &END" << endl;
89125dc7 1002 out << name.Data() << endl;
1003 for (i=0; i<nel; i++) {
1004 elem = mix->GetElement(i);
1005 line = elem->GetName();
1006 if (line.Length()==1) line.Append(" ");
1007 out << line.Data() << " ";
1008 }
1009 out << endl;
95161cd7 1010
1011 delete [] zt;
1012 delete [] wt;
89125dc7 1013 }
95161cd7 1014
6d184c54 1015 Double_t ue = 3000000.; // [MeV]
1016 Double_t up = 3000000.; // [MeV]
1017 Double_t ae = -1.;
1018 Double_t ap = -1.;
1019
1020
1021 TObjArray* cutList = ((TFluka*) gMC)->GetListOfUserConfigs();
1022 TIter next(cutList);
1023 TFlukaConfigOption* proc;
1024
1025 while((proc = (TFlukaConfigOption*)next()))
1026 {
1027 if (proc->Medium() == mat->GetIndex()) {
4aba9d66 1028 ap = proc->Cut(kCUTGAM);
1029 ae = proc->Cut(kCUTELE);
1030 if (ap == -1.) ap = TFlukaConfigOption::DefaultCut(kCUTGAM);
1031 if (ae == -1.) ae = TFlukaConfigOption::DefaultCut(kCUTELE);
1032 break;
6d184c54 1033 }
1034 }
1035
1036 if (ap == -1.) ap = TFlukaConfigOption::DefaultCut(kCUTGAM);
1037 if (ae == -1.) ae = TFlukaConfigOption::DefaultCut(kCUTELE);
1038
1039 ap *= 1000.; // [MeV]
1040 ae = (ae + 0.00051099906) * 1000.; // [MeV]
1041
89125dc7 1042 out << "ENER" << endl;
6d184c54 1043 out << " $INP AE=" << ae << ", UE=" << ue <<", AP=" << ap << ", UP=" << up << " $END" << endl;
89125dc7 1044 out << "PWLF" << endl;
1045 out << " $INP NALE=300, NALG=400, NALR=100 $END" << endl;
1046 out << "DECK" << endl;
1047 out << " $INP $END" << endl;
1048 out << "TEST" << endl;
1049 out << " $INP $END" << endl;
1050 out.close();
8495a208 1051}
1052
95161cd7 1053Double_t * TFlukaMCGeometry::GetISSB(Double_t rho, Int_t nElem, Double_t *zelem, Double_t *welem ) const
1054{
1055 // Read the density effect parameters
1056 // from R.M. Sternheimer et al. Atomic Data
1057 // and Nuclear Data Tables, Vol. 30 No. 2
1058 //
1059 // return the parameters if the element/mixture match with one of the list
1060 // otherwise returns the parameters set to 0
1061
1062 struct sternheimerData {
1063 TString longname; // element/mixture name
1064 Int_t nelems; // number of constituents N
1065 Int_t Z[20]; //[nelems] Z
1066 Double_t wt[20]; //[nelems] weight fraction
1067 Double_t density; // g/cm3
1068 Double_t iev; // Average Ion potential (eV)
6d184c54 1069 // **** Sternheimer parameters ****
95161cd7 1070 Double_t cbar; // CBAR
1071 Double_t x0; // X0
1072 Double_t x1; // X1
1073 Double_t afact; // AFACT
1074 Double_t sk; // SK
1075 Double_t delta0; // DELTA0
4aba9d66 1076
1077 sternheimerData():
1078 longname(""), nelems(0), density(0), iev(0), cbar(0),
1079 x0(0), x1(0), afact(0), sk(0), delta0(0) {}
95161cd7 1080 };
1081
1082 TString shortname;
1083 TString formula;
1084 Int_t num;
1085 char state;
1086
1087 static Double_t parameters[7];
1088 memset( parameters, 0, sizeof(Double_t) );
1089
1090 static sternheimerData sternDataArray[300];
1091 static Bool_t isFileRead = kFALSE;
1092
1093 // Read the data file if is needed
1094 if( isFileRead == kFALSE ) {
1095 TString sSternheimerInp = getenv("ALICE_ROOT");
1096 sSternheimerInp +="/TFluka/input/Sternheimer.data";
1097
1098 ifstream in(sSternheimerInp);
1099 char line[100];
1100 in.getline(line, 100);
1101 in.getline(line, 100);
1102 in.getline(line, 100);
1103 in.getline(line, 100);
1104 in.getline(line, 100);
1105 in.getline(line, 100);
1106
1107
1108 Int_t is = 0;
1109 while( !in.eof() ) {
1110 in >> shortname >> num >> sternDataArray[is].nelems
1111 >> sternDataArray[is].longname >> formula >> state;
1112 if( in.eof() ) break;
1113 for(int i=0; i<sternDataArray[is].nelems; i++) {
1114 in >> sternDataArray[is].Z[i] >> sternDataArray[is].wt[i];
1115 }
1116 in >> sternDataArray[is].density;
1117 in >> sternDataArray[is].iev;
1118 in >> sternDataArray[is].cbar;
1119 in >> sternDataArray[is].x0;
1120 in >> sternDataArray[is].x1;
1121 in >> sternDataArray[is].afact;
1122 in >> sternDataArray[is].sk;
1123 if( sternDataArray[is].nelems == 1 ) in >> sternDataArray[is].delta0;
1124 is++;
1125 }
1126 isFileRead = kTRUE;
1127 in.close();
1128 }
1129
1130 Int_t is = 0;
1131 while( is < 280 ) {
1132
1133 // check for elements
1134 if( sternDataArray[is].nelems == 1 && nElem == 1
1135 && sternDataArray[is].Z[0] == Int_t(*zelem)
1136 && TMath::Abs( (sternDataArray[is].density - rho)/sternDataArray[is].density ) < 0.1 ) {
1137 cout << sternDataArray[is].longname << " #elems:" << sternDataArray[is].nelems << " Rho:"
1138 << sternDataArray[is].density << endl;
1139 cout << sternDataArray[is].iev << " "
1140 << sternDataArray[is].cbar << " "
1141 << sternDataArray[is].x0 << " "
1142 << sternDataArray[is].x1 << " "
1143 << sternDataArray[is].afact << " "
1144 << sternDataArray[is].sk << " "
1145 << sternDataArray[is].delta0 << endl;
1146
1147 parameters[0] = sternDataArray[is].iev;
1148 parameters[1] = sternDataArray[is].cbar;
1149 parameters[2] = sternDataArray[is].x0;
1150 parameters[3] = sternDataArray[is].x1;
1151 parameters[4] = sternDataArray[is].afact;
1152 parameters[5] = sternDataArray[is].sk;
1153 parameters[6] = sternDataArray[is].delta0;
1154 return parameters;
1155 }
1156
1157 // check for mixture
1158 int nmatch = 0;
1159 if( sternDataArray[is].nelems > 1 && sternDataArray[is].nelems == nElem ) {
1160 for(int j=0; j<sternDataArray[is].nelems; j++) {
1161 if( sternDataArray[is].Z[j] == Int_t(zelem[j]) &&
1162 TMath::Abs( (sternDataArray[is].wt[j] - welem[j])/sternDataArray[is].wt[j] ) < 0.1 )
1163 nmatch++;
1164 }
1165 }
1166
1167 if( sternDataArray[is].nelems > 1 &&
1168 TMath::Abs( (sternDataArray[is].density - rho)/sternDataArray[is].density ) < 0.1
1169 && nmatch == sternDataArray[is].nelems ) {
1170 cout << sternDataArray[is].longname << " #elem:" << sternDataArray[is].nelems << " Rho:"
1171 << sternDataArray[is].density << endl;
1172 cout << sternDataArray[is].iev << " "
1173 << sternDataArray[is].cbar << " "
1174 << sternDataArray[is].x0 << " "
1175 << sternDataArray[is].x1 << " "
1176 << sternDataArray[is].afact << " "
1177 << sternDataArray[is].sk << " "
1178 << sternDataArray[is].delta0 << endl;
1179
1180 parameters[0] = sternDataArray[is].iev;
1181 parameters[1] = sternDataArray[is].cbar;
1182 parameters[2] = sternDataArray[is].x0;
1183 parameters[3] = sternDataArray[is].x1;
1184 parameters[4] = sternDataArray[is].afact;
1185 parameters[5] = sternDataArray[is].sk;
1186 parameters[6] = 0;
1187 return parameters;
1188 }
1189 is++;
1190 }
1191 return parameters;
1192}
1193
8495a208 1194//_____________________________________________________________________________
1195void TFlukaMCGeometry::PrintHeader(ofstream &out, const char *text) const
1196{
1197// Print a FLUKA header.
1198 out << "*\n" << "*\n" << "*\n";
1199 out << "********************* " << text << " *********************\n"
1200 << "*\n";
1201 out << "*...+....1....+....2....+....3....+....4....+....5....+....6....+....7..."
1202 << endl;
1203 out << "*" << endl;
1204}
1205
1206//_____________________________________________________________________________
1207Int_t TFlukaMCGeometry::RegionId() const
1208{
1209// Returns current region id <-> TGeo node id
1210 if (gGeoManager->IsOutside()) return 0;
1211 return gGeoManager->GetCurrentNode()->GetUniqueID();
1212}
89125dc7 1213
1214//_____________________________________________________________________________
1215Int_t TFlukaMCGeometry::GetElementIndex(Int_t z) const
1216{
6f5667d1 1217// Get index of a material having a given Z element.
89125dc7 1218 TIter next(fMatList);
1219 TGeoMaterial *mat;
1220 Int_t index = 0;
1221 while ((mat=(TGeoMaterial*)next())) {
1222 if (mat->IsMixture()) continue;
1223 if (mat->GetElement()->Z() == z) return mat->GetIndex();
1224 }
1225 return index;
1226}
1227
05265ca9 1228//_____________________________________________________________________________
aae6bcdd 1229void TFlukaMCGeometry::SetMreg(Int_t mreg, Int_t lttc)
05265ca9 1230{
1231// Update if needed next history;
aae6bcdd 1232// if (gFluka->GetDummyBoundary()==2) {
1233// gGeoManager->CdNode(fNextLattice-1);
1234// return;
1235// }
1236 if (lttc == TFlukaMCGeometry::kLttcOutside) {
1237 fCurrentRegion = NofVolumes()+2;
1238 fCurrentLattice = lttc;
1239 gGeoManager->CdTop();
1240 gGeoManager->SetOutside(kTRUE);
1241 }
1242 if (lttc == TFlukaMCGeometry::kLttcVirtual) return;
1243 if (lttc <=0) {
1244 Error("TFlukaMCGeometry::SetMreg","Invalide lattice %i",lttc);
05265ca9 1245 return;
aae6bcdd 1246 }
1247 fCurrentRegion = mreg;
1248 fCurrentLattice = lttc;
1249
1250 Int_t crtlttc = gGeoManager->GetCurrentNodeId()+1;
1251 if (crtlttc == lttc) return;
1252 gGeoManager->CdNode(lttc-1);
9968e86c 1253 while (gGeoManager->GetCurrentVolume()->IsAssembly()) gGeoManager->CdUp();
05265ca9 1254}
1255
1256//_____________________________________________________________________________
d11bf3ca 1257void TFlukaMCGeometry::SetCurrentRegion(Int_t mreg, Int_t latt)
1258{
1259// Set index/history for next entered region
1260 fCurrentRegion = mreg;
1261 fCurrentLattice = latt;
1262}
1263
1264//_____________________________________________________________________________
05265ca9 1265void TFlukaMCGeometry::SetNextRegion(Int_t mreg, Int_t latt)
1266{
1267// Set index/history for next entered region
1268 fNextRegion = mreg;
1269 fNextLattice = latt;
1270}
8495a208 1271
1272//_____________________________________________________________________________
1273void TFlukaMCGeometry::ToFlukaString(TString &str) const
1274{
1275// ToFlukaString converts an string to something usefull in FLUKA:
1276// * Capital letters
1277// * Only 8 letters
1278// * Replace ' ' by '_'
1279 if (str.Length()<8) {
1280 str += " ";
1281 }
1282 str.Remove(8);
1283 Int_t ilast;
1284 for (ilast=7; ilast>0; ilast--) if (str(ilast)!=' ') break;
1285 str.ToUpper();
1286 for (Int_t pos=0; pos<ilast; pos++)
1287 if (str(pos)==' ') str.Replace(pos,1,"_",1);
1288 return;
1289}
89125dc7 1290
1291//_____________________________________________________________________________
1292void TFlukaMCGeometry::FlukaMatName(TString &str) const
1293{
81f1d030 1294// Strip the detector name
1295
1296 TObjArray * tokens = str.Tokenize("_");
1297 Int_t ntok = tokens->GetEntries();
1298 if (ntok > 0) {
4aba9d66 1299 TString head = ((TObjString*) tokens->At(0))->GetString();
1300 Int_t nhead = head.Length();
1301 str = str.Remove(0, nhead + 1);
81f1d030 1302 }
1303 tokens->Clear();
1304 delete tokens;
1305
6f5667d1 1306// Convert a name to upper case 8 chars.
89125dc7 1307 ToFlukaString(str);
1308 Int_t ilast;
1309 for (ilast=7; ilast>0; ilast--) if (str(ilast)!=' ') break;
1310 if (ilast>5) ilast = 5;
1311 char number[3];
1312 TIter next(fMatNames);
1313 TObjString *objstr;
1314 TString matname;
1315 Int_t index = 0;
1316 while ((objstr=(TObjString*)next())) {
1317 matname = objstr->GetString();
1318 if (matname == str) {
1319 index++;
1320 if (index<10) {
1321 number[0] = '0';
1322 sprintf(&number[1], "%d", index);
1323 } else if (index<100) {
1324 sprintf(number, "%d", index);
1325 } else {
1326 Error("FlukaMatName", "Too many materials %s", str.Data());
1327 return;
1328 }
1329 str.Replace(ilast+1, 2, number);
1330 str.Remove(8);
1331 }
1332 }
1333}
4aba9d66 1334
8495a208 1335//______________________________________________________________________________
1336void TFlukaMCGeometry::Vname(const char *name, char *vname) const
1337{
1338 //
1339 // convert name to upper case. Make vname at least 4 chars
1340 //
1341 Int_t l = strlen(name);
1342 Int_t i;
1343 l = l < 4 ? l : 4;
1344 for (i=0;i<l;i++) vname[i] = toupper(name[i]);
1345 for (i=l;i<4;i++) vname[i] = ' ';
1346 vname[4] = 0;
1347}
1348
1349
4aba9d66 1350//______________________________________________________________________________
1351Int_t TFlukaMCGeometry::GetNstep()
1352{
1353 // return gNstep for debug propose
1354 return gNstep;
1355}
81f1d030 1356
8495a208 1357
9968e86c 1358Int_t TFlukaMCGeometry::GetPredefinedMaterialId(Int_t z) const
1359{
1360// Get predifined material id from Z if present in list
1361 const Int_t kMax = 25;
1362
1363 static Int_t idFluka[kMax] =
1364 {-1, -1, 1, 2, 4, 6, 7, 8, 12, 13,
1365 26, 29, 47, 14, 79, 80, 82, 73, 11, 18,
1366 20, 50, 74, 22, 28};
1367
1368 Int_t id = -1;
1369
1370 for (Int_t i = 0; i < kMax; i++)
1371 {
1372 if (z == idFluka[i]) {
1373 id = i + 1;
1374 break;
1375 }
1376
1377 }
1378
1379 return id;
1380}
1381
1382// FLUKA GEOMETRY WRAPPERS - to replace FLUGG wrappers
1383//
8495a208 1384//_____________________________________________________________________________
8495a208 1385Int_t idnrwr(const Int_t & /*nreg*/, const Int_t & /*mlat*/)
1386{
1387// from FLUGG:
1388// Wrapper for setting DNEAR option on fluka side. Must return 0
1389// if user doesn't want Fluka to use DNEAR to compute the
1390// step (the same effect is obtained with the GLOBAL (WHAT(3)=-1)
1391// card in fluka input), returns 1 if user wants Fluka always to
1392// use DNEAR (in this case, be sure that GEANT4 DNEAR is unique,
1393// coming from all directions!!!)
89125dc7 1394 if (gMCGeom->IsDebugging()) printf("========== Dummy IDNRWR\n");
8495a208 1395 return 0;
1396}
1397
1398//_____________________________________________________________________________
8495a208 1399void g1wr(Double_t &pSx, Double_t &pSy, Double_t &pSz,
d11bf3ca 1400 Double_t *pV, Int_t &oldReg , const Int_t &oldLttc, Double_t &propStep,
cefd9d9b 1401 Int_t &/*nascFlag*/, Double_t &retStep, Int_t &newReg,
4aba9d66 1402 Double_t &saf, Int_t &newLttc, Int_t &lttcFlag,
8495a208 1403 Double_t *sLt, Int_t *jrLt)
d11bf3ca 1404
8495a208 1405{
d11bf3ca 1406 // Initialize FLUKa point and direction;
4b1c3976 1407 static Int_t ierr = 0;
6f5667d1 1408 gNstep++;
4b1c3976 1409// gMCGeom->SetDebugMode(kTRUE);
4aba9d66 1410
784c3bdd 1411 NORLAT.xn[0] = pSx;
1412 NORLAT.xn[1] = pSy;
1413 NORLAT.xn[2] = pSz;
8495a208 1414
a940d455 1415 Int_t olttc = oldLttc;
81f1d030 1416 if (oldLttc<=0) {
a940d455 1417 gGeoManager->FindNode(pSx,pSy,pSz);
1418 olttc = gGeoManager->GetCurrentNodeId()+1;
4aba9d66 1419 oldReg = gGeoManager->GetCurrentVolume()->GetNumber();
1420 }
1421
4b1c3976 1422 if (gMCGeom->IsDebugging()) {
1423 cout << "g1wr gNstep=" << gNstep << " oldReg="<< oldReg <<" olttc="<< olttc
1424 << " track=" << TRACKR.ispusr[mkbmx2-1] << endl;
1425 cout << " point: (" << pSx << ", " << pSy << ", " << pSz << ") dir: ("
1426 << pV[0] << ", " << pV[1] << ", " << pV[2] << ")" << endl;
1427 }
1428
1429
1430 Int_t ccreg=0,cclat=0;
81f1d030 1431 gMCGeom->GetCurrentRegion(ccreg,cclat);
4aba9d66 1432 Bool_t crossed = (ccreg==oldReg && cclat==olttc)?kFALSE:kTRUE;
1433
a940d455 1434 gMCGeom->SetCurrentRegion(oldReg, olttc);
d11bf3ca 1435 // Initialize default return values
1436 lttcFlag = 0;
a940d455 1437 jrLt[lttcFlag] = olttc;
d11bf3ca 1438 sLt[lttcFlag] = propStep;
1439 jrLt[lttcFlag+1] = -1;
1440 sLt[lttcFlag+1] = 0.;
1441 newReg = oldReg;
a940d455 1442 newLttc = olttc;
aae6bcdd 1443 Bool_t crossedDummy = (oldReg == gFluka->GetDummyRegion())?kTRUE:kFALSE;
d11bf3ca 1444 Int_t curLttc, curReg;
aae6bcdd 1445 if (crossedDummy) {
1446 // FLUKA crossed the dummy boundary - update new region/history
1447 retStep = TGeoShape::Tolerance();
1448 saf = 0.;
1449 gMCGeom->GetNextRegion(newReg, newLttc);
1450 gMCGeom->SetMreg(newReg, newLttc);
1451 sLt[lttcFlag] = TGeoShape::Tolerance(); // null step in current region
1452 lttcFlag++;
1453 jrLt[lttcFlag] = newLttc;
1454 sLt[lttcFlag] = TGeoShape::Tolerance(); // null step in next region
1455 jrLt[lttcFlag+1] = -1;
1456 sLt[lttcFlag+1] = 0.; // null step in next region;
4aba9d66 1457 if (gMCGeom->IsDebugging()) printf("=> crossed dummy!! newReg=%i newLttc=%i\n", newReg, newLttc);
aae6bcdd 1458 return;
4aba9d66 1459 }
1460
d11bf3ca 1461 // Reset outside flag
cefd9d9b 1462 gGeoManager->SetOutside(kFALSE);
4aba9d66 1463
d11bf3ca 1464 curLttc = gGeoManager->GetCurrentNodeId()+1;
1465 curReg = gGeoManager->GetCurrentVolume()->GetNumber();
a940d455 1466 if (olttc != curLttc) {
d11bf3ca 1467 // FLUKA crossed the boundary : we trust that the given point is really there,
1468 // so we just update TGeo state
a940d455 1469 gGeoManager->CdNode(olttc-1);
2573ac89 1470 curLttc = gGeoManager->GetCurrentNodeId()+1;
d11bf3ca 1471 curReg = gGeoManager->GetCurrentVolume()->GetNumber();
4aba9d66 1472 }
784c3bdd 1473 // Now the current TGeo state reflects the FLUKA state
4aba9d66 1474
aae6bcdd 1475 gGeoManager->SetCurrentPoint(pSx, pSy, pSz);
cefd9d9b 1476 gGeoManager->SetCurrentDirection(pV);
4b1c3976 1477 Double_t pt[3], local[3], ldir[3];
1478 Int_t pid = TRACKR.jtrack;
1479 pt[0] = pSx;
1480 pt[1] = pSy;
1481 pt[2] = pSz;
1482 gGeoManager->MasterToLocal(pt,local);
1483 gGeoManager->MasterToLocalVect(pV,ldir);
1484/*
1485 Bool_t valid = gGeoManager->GetCurrentVolume()->Contains(local);
1486 if (!valid) {
1487 printf("location not valid in %s pid=%i\n", gGeoManager->GetPath(),pid);
1488 printf("local:(%f, %f, %f) ldir:(%f, %f, %f)\n", local[0],local[1],local[2],ldir[0],ldir[1],ldir[2]);
1489// gGeoManager->FindNode();
1490// printf(" -> actual location: %s\n", gGeoManager->GetPath());
1491 }
1492*/
1493 Double_t pstep = propStep;
1494 Double_t snext = propStep;
1495 const Double_t epsil = 0.9999999 * TGeoShape::Tolerance();
1496 // This should never happen !!!
1497 if (pstep<TGeoShape::Tolerance()) {
1498 printf("Proposed step is 0 !!!\n");
1499 pstep = 2.*TGeoShape::Tolerance();
1500 }
aae6bcdd 1501 if (crossed) {
4b1c3976 1502 gGeoManager->FindNextBoundaryAndStep(pstep);
1503 snext = gGeoManager->GetStep();
1504 saf = 0.;
1505 if (snext <= TGeoShape::Tolerance()) {
1506// printf("FLUKA crossed boundary but snext=%g\n", snext);
1507 ierr++;
1508 snext = epsil;
1509 } else {
1510 snext += TGeoShape::Tolerance();
1511 ierr = 0;
1512 }
aae6bcdd 1513 } else {
4b1c3976 1514 gGeoManager->FindNextBoundaryAndStep(pstep, kTRUE);
1515 snext = gGeoManager->GetStep();
aae6bcdd 1516 saf = gGeoManager->GetSafeDistance();
4b1c3976 1517 if (snext <= TGeoShape::Tolerance()) {
1518// printf("FLUKA put particle on bondary without crossing\n");
1519 ierr++;
1520 snext = epsil;
1521 saf = 0.;
1522 } else {
1523 snext += TGeoShape::Tolerance();
1524 ierr = 0;
1525 }
1526 if (saf<0) saf=0.;
aae6bcdd 1527 saf -= saf*3.0e-09;
4aba9d66 1528 }
4b1c3976 1529// if (ierr>1) {
1530// printf("%d snext=%g\n", ierr, snext);
1531// }
1532 if (ierr == 10) {
1533// printf("Too many null steps - sending error code -33...\n");
1534 newReg = -33; // Error code
1535 ierr = 0;
1536 return;
1537 }
1538
81f1d030 1539 PAREM.dist = snext;
784c3bdd 1540 NORLAT.distn = snext;
1541 NORLAT.xn[0] += snext*pV[0];
1542 NORLAT.xn[1] += snext*pV[1];
1543 NORLAT.xn[2] += snext*pV[2];
aae6bcdd 1544 if (!gGeoManager->IsOnBoundary()) {
d11bf3ca 1545 // Next boundary further than proposed step, which is approved
4aba9d66 1546 if (saf>propStep) saf = propStep;
d11bf3ca 1547 retStep = propStep;
1548 sLt[lttcFlag] = propStep;
1549 return;
1550 }
aae6bcdd 1551 if (saf>snext) saf = snext; // Safety should be less than the proposed step if a boundary will be crossed
cefd9d9b 1552 gGeoManager->SetCurrentPoint(pSx,pSy,pSz);
d11bf3ca 1553 newLttc = (gGeoManager->IsOutside())?(TFlukaMCGeometry::kLttcOutside):gGeoManager->GetCurrentNodeId()+1;
aae6bcdd 1554 newReg = (gGeoManager->IsOutside())?(gMCGeom->NofVolumes()+2):gGeoManager->GetCurrentVolume()->GetNumber();
4aba9d66 1555 if (gMCGeom->IsDebugging()) printf("=> newReg=%i newLttc=%i\n", newReg, newLttc);
d11bf3ca 1556
1557 // We really crossed the boundary, but is it the same region ?
89125dc7 1558 gMCGeom->SetNextRegion(newReg, newLttc);
aae6bcdd 1559
4aba9d66 1560 if ( ((newReg==oldReg && newLttc!=olttc) || (oldReg!=newReg && olttc==newLttc) ) && pid!=-1) {
d11bf3ca 1561 // Virtual boundary between replicants
aae6bcdd 1562 newReg = gFluka->GetDummyRegion();
d11bf3ca 1563 newLttc = TFlukaMCGeometry::kLttcVirtual;
4aba9d66 1564 if (gMCGeom->IsDebugging()) printf("=> virtual boundary!! newReg=%i newLttc=%i\n", newReg, newLttc);
1565 }
1566
d11bf3ca 1567 retStep = snext;
1568 sLt[lttcFlag] = snext;
1569 lttcFlag++;
1570 jrLt[lttcFlag] = newLttc;
1571 sLt[lttcFlag] = snext;
1572 jrLt[lttcFlag+1] = -1;
4aba9d66 1573
1574 sLt[lttcFlag+1] = 0.;
87f4b393 1575 gGeoManager->SetOutside(kFALSE);
1576 gGeoManager->CdNode(olttc-1);
89125dc7 1577 if (gMCGeom->IsDebugging()) {
d11bf3ca 1578 printf("=> snext=%g safe=%g\n", snext, saf);
1579 for (Int_t i=0; i<lttcFlag+1; i++) printf(" jrLt[%i]=%i sLt[%i]=%g\n", i,jrLt[i],i,sLt[i]);
4aba9d66 1580 }
8495a208 1581}
1582
efde9b4d 1583//_____________________________________________________________________________
1584void g1rtwr()
1585{
4aba9d66 1586
89125dc7 1587 if (gMCGeom->IsDebugging()) printf("========== Dummy G1RTWR\n");
efde9b4d 1588}
1589
1590//_____________________________________________________________________________
4aba9d66 1591void conhwr(Int_t & intHist, Int_t & incrCount)
efde9b4d 1592{
4aba9d66 1593 if (gMCGeom->IsDebugging()) printf("========== Dummy CONHWR intHist=%d incrCount=%d currentNodeId=%d\n",
1594 intHist, incrCount, gGeoManager->GetCurrentNodeId()+1 );
1595// if( incrCount != -1 ) {
1596// if (intHist==0) gGeoManager->CdTop();
1597// else gGeoManager->CdNode(intHist-1);
1598// }
1599// intHist = gGeoManager->GetCurrentNodeId()+1;
efde9b4d 1600}
1601
1602//_____________________________________________________________________________
2573ac89 1603void inihwr(Int_t &intHist)
efde9b4d 1604{
4aba9d66 1605 if (gMCGeom->IsDebugging())
1606 printf("========== Inside INIHWR -> reinitializing history: %i \n", intHist);
2573ac89 1607 if (gGeoManager->IsOutside()) gGeoManager->CdTop();
4aba9d66 1608 if (intHist<0) {
2573ac89 1609// printf("=== wrong history number\n");
1610 return;
1611 }
1612 if (intHist==0) gGeoManager->CdTop();
1613 else gGeoManager->CdNode(intHist-1);
89125dc7 1614 if (gMCGeom->IsDebugging()) {
2bc4c610 1615 printf(" --- current path: %s\n", gGeoManager->GetPath());
1616 printf("<= INIHWR\n");
4aba9d66 1617 }
efde9b4d 1618}
1619
1620//_____________________________________________________________________________
1621void jomiwr(const Int_t & /*nge*/, const Int_t & /*lin*/, const Int_t & /*lou*/,
1622 Int_t &flukaReg)
1623{
1624// Geometry initialization wrapper called by FLUKAM. Provides to FLUKA the
1625// number of regions (volumes in TGeo)
1626 // build application geometry
89125dc7 1627 if (gMCGeom->IsDebugging()) printf("========== Inside JOMIWR\n");
aae6bcdd 1628 flukaReg = gGeoManager->GetListOfUVolumes()->GetEntriesFast()+1;
89125dc7 1629 if (gMCGeom->IsDebugging()) printf("<= JOMIWR: last region=%i\n", flukaReg);
efde9b4d 1630}
1631
1632//_____________________________________________________________________________
1633void lkdbwr(Double_t &pSx, Double_t &pSy, Double_t &pSz,
784c3bdd 1634 Double_t *pV, const Int_t &oldReg, const Int_t &oldLttc,
aae6bcdd 1635 Int_t &flagErr, Int_t &newReg, Int_t &newLttc)
efde9b4d 1636{
89125dc7 1637 if (gMCGeom->IsDebugging()) {
2bc4c610 1638 printf("========== Inside LKDBWR (%f, %f, %f)\n",pSx, pSy, pSz);
784c3bdd 1639 printf(" in: pV=(%f, %f, %f)\n", pV[0], pV[1], pV[2]);
2bc4c610 1640 printf(" in: oldReg=%i oldLttc=%i\n", oldReg, oldLttc);
1641 }
aae6bcdd 1642 lkwr(pSx,pSy,pSz,pV,oldReg,oldLttc,flagErr,newReg,newLttc);
efde9b4d 1643}
1644
1645//_____________________________________________________________________________
1646void lkfxwr(Double_t &pSx, Double_t &pSy, Double_t &pSz,
784c3bdd 1647 Double_t *pV, const Int_t &oldReg, const Int_t &oldLttc,
aae6bcdd 1648 Int_t &flagErr, Int_t &newReg, Int_t &newLttc)
efde9b4d 1649{
89125dc7 1650 if (gMCGeom->IsDebugging()) {
2bc4c610 1651 printf("========== Inside LKFXWR (%f, %f, %f)\n",pSx, pSy, pSz);
784c3bdd 1652 printf(" in: pV=(%f, %f, %f)\n", pV[0], pV[1], pV[2]);
2bc4c610 1653 printf(" in: oldReg=%i oldLttc=%i\n", oldReg, oldLttc);
1654 }
aae6bcdd 1655 lkwr(pSx,pSy,pSz,pV,oldReg,oldLttc,flagErr,newReg,newLttc);
efde9b4d 1656}
1657
1658//_____________________________________________________________________________
1659void lkmgwr(Double_t &pSx, Double_t &pSy, Double_t &pSz,
784c3bdd 1660 Double_t *pV, const Int_t &oldReg, const Int_t &oldLttc,
4aba9d66 1661 Int_t &flagErr, Int_t &newReg, Int_t &newLttc)
efde9b4d 1662{
89125dc7 1663 if (gMCGeom->IsDebugging()) {
2bc4c610 1664 printf("========== Inside LKMGWR (%f, %f, %f)\n",pSx, pSy, pSz);
784c3bdd 1665 printf(" in: pV=(%f, %f, %f)\n", pV[0], pV[1], pV[2]);
2bc4c610 1666 printf(" in: oldReg=%i oldLttc=%i\n", oldReg, oldLttc);
1667 }
aae6bcdd 1668 lkwr(pSx,pSy,pSz,pV,oldReg,oldLttc,flagErr,newReg,newLttc);
efde9b4d 1669}
1670
1671//_____________________________________________________________________________
1672void lkwr(Double_t &pSx, Double_t &pSy, Double_t &pSz,
784c3bdd 1673 Double_t *pV, const Int_t &oldReg, const Int_t &oldLttc,
4aba9d66 1674 Int_t &flagErr, Int_t &newReg, Int_t &newLttc)
efde9b4d 1675{
89125dc7 1676 if (gMCGeom->IsDebugging()) {
2bc4c610 1677 printf("========== Inside LKWR (%f, %f, %f)\n",pSx, pSy, pSz);
784c3bdd 1678 printf(" in: pV=(%f, %f, %f)\n", pV[0], pV[1], pV[2]);
2bc4c610 1679 printf(" in: oldReg=%i oldLttc=%i\n", oldReg, oldLttc);
1680 }
aae6bcdd 1681 flagErr = 0;
4b1c3976 1682 Double_t epsil = 1000.*TGeoShape::Tolerance();
1683 TGeoNode *node = gGeoManager->FindNode(pSx+epsil*pV[0], pSy+epsil*pV[1], pSz+epsil*pV[2]);
efde9b4d 1684 if (gGeoManager->IsOutside()) {
aae6bcdd 1685 newReg = gMCGeom->NofVolumes()+2;
1686 newLttc = TFlukaMCGeometry::kLttcOutside;
87f4b393 1687 gGeoManager->SetOutside(kFALSE);
1688 if (oldLttc>0 && oldLttc<newLttc) gGeoManager->CdNode(oldLttc-1);
05265ca9 1689 return;
efde9b4d 1690 }
87f4b393 1691 gGeoManager->SetOutside(kFALSE);
efde9b4d 1692 newReg = node->GetVolume()->GetNumber();
d11bf3ca 1693 newLttc = gGeoManager->GetCurrentNodeId()+1;
aae6bcdd 1694 if (oldLttc==TFlukaMCGeometry::kLttcOutside || oldLttc==0) return;
1695
1696 Int_t dummy = gFluka->GetDummyRegion();
1697 if (oldReg==dummy) {
1698 Int_t newreg1, newlttc1;
1699 gMCGeom->GetNextRegion(newreg1, newlttc1);
1700 if (newreg1==newReg && newlttc1==newLttc) {
1701 newReg = dummy;
1702 newLttc = TFlukaMCGeometry::kLttcVirtual;
4aba9d66 1703 flagErr = newReg;
1704 if (gMCGeom->IsDebugging()) printf(" virtual boundary (oldReg==dummy) !! newReg=%i newLttc=%i\n", newReg, newLttc);
1705 }
aae6bcdd 1706 return;
4aba9d66 1707 }
aae6bcdd 1708
1709 if (oldReg==newReg && oldLttc!=newLttc) {
4aba9d66 1710 newReg = dummy;
aae6bcdd 1711 newLttc = TFlukaMCGeometry::kLttcVirtual;
4aba9d66 1712 if (gMCGeom->IsDebugging()) printf(" virtual boundary!! newReg=%i newLttc=%i\n", newReg, newLttc);
1713 }
1714
1715 if( oldReg!=newReg && oldLttc==newLttc ) {
1716 // this should not happen!! ??? Ernesto
4b1c3976 1717// cout << " lkwr oldReg!=newReg ("<< oldReg <<"!="<< newReg
1718// << ") && oldLttc==newLttc ("<< newLttc <<") !!!!!!!!!!!!!!!!!" << endl;
4aba9d66 1719 newReg = dummy;
1720 newLttc = TFlukaMCGeometry::kLttcVirtual;
1721 flagErr = newReg;
1722 }
1723
89125dc7 1724 if (gMCGeom->IsDebugging()) {
aae6bcdd 1725 printf(" LKWR: newReg=%i newLttc=%i\n", newReg, newLttc);
2bc4c610 1726 }
efde9b4d 1727}
1728
1729//_____________________________________________________________________________
2573ac89 1730void nrmlwr(Double_t &pSx, Double_t &pSy, Double_t &pSz,
1731 Double_t &pVx, Double_t &pVy, Double_t &pVz,
4aba9d66 1732 Double_t *norml, const Int_t &oldReg,
1733 const Int_t &newReg, Int_t &flagErr)
efde9b4d 1734{
89125dc7 1735 if (gMCGeom->IsDebugging()) {
2bc4c610 1736 printf("========== Inside NRMLWR (%g, %g, %g, %g, %g, %g)\n", pSx,pSy,pSz,pVx,pVy,pVz);
4aba9d66 1737 printf(" (%g, %g, %g)\n", NORLAT.xn[0], NORLAT.xn[1], NORLAT.xn[2]);
2bc4c610 1738 printf(" oldReg=%i, newReg=%i\n", oldReg,newReg);
1739 }
784c3bdd 1740 gGeoManager->SetCurrentPoint(NORLAT.xn[0], NORLAT.xn[1], NORLAT.xn[2]);
aae6bcdd 1741 gGeoManager->SetCurrentDirection(pVx, pVy, pVz);
2573ac89 1742 Double_t *dnorm = gGeoManager->FindNormalFast();
1743 flagErr = 0;
1744 if (!dnorm) {
1745 printf(" ERROR: Cannot compute fast normal\n");
1746 flagErr = 1;
1747 norml[0] = -pVx;
1748 norml[1] = -pVy;
1749 norml[2] = -pVz;
aae6bcdd 1750 } else {
4aba9d66 1751 norml[0] = -dnorm[0];
1752 norml[1] = -dnorm[1];
1753 norml[2] = -dnorm[2];
aae6bcdd 1754 }
4aba9d66 1755
89125dc7 1756 if (gMCGeom->IsDebugging()) {
784c3bdd 1757 printf(" normal to boundary: (%g, %g, %g)\n", norml[0], norml[1], norml[2]);
2bc4c610 1758 printf("<= NRMLWR\n");
4aba9d66 1759 }
1760
efde9b4d 1761}
1762
1763//_____________________________________________________________________________
1764void rgrpwr(const Int_t & /*flukaReg*/, const Int_t & /*ptrLttc*/, Int_t & /*g4Reg*/,
1765 Int_t * /*indMother*/, Int_t * /*repMother*/, Int_t & /*depthFluka*/)
1766{
89125dc7 1767 if (gMCGeom->IsDebugging()) printf("=> Dummy RGRPWR\n");
efde9b4d 1768}
1769
1770//_____________________________________________________________________________
1771Int_t isvhwr(const Int_t &check, const Int_t & intHist)
1772{
1773// from FLUGG:
1774// Wrapper for saving current navigation history (fCheck=default)
1775// and returning its pointer. If fCheck=-1 copy of history pointed
1776// by intHist is made in NavHistWithCount object, and its pointer
1777// is returned. fCheck=1 and fCheck=2 cases are only in debugging
1778// version: an array is created by means of FGeometryInit functions
1779// (but could be a static int * ptrArray = new int[10000] with
1780// file scope as well) that stores a flag for deleted/undeleted
1781// histories and at the end of event is checked to verify that
1782// all saved history objects have been deleted.
1783
1784// For TGeo, just return the current node ID. No copy need to be made.
1785
4aba9d66 1786 if (gMCGeom->IsDebugging()) printf("=> Inside ISVHWR check=%d intHist=%d\n", check, intHist);
efde9b4d 1787 if (check<0) return intHist;
2573ac89 1788 Int_t histInt = gGeoManager->GetCurrentNodeId()+1;
89125dc7 1789 if (gMCGeom->IsDebugging()) printf("<= ISVHWR: history is: %i in: %s\n", histInt, gGeoManager->GetPath());
efde9b4d 1790 return histInt;
1791}
1792
1793
8495a208 1794
1795