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