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