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