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