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