]> git.uio.no Git - u/mrichter/AliRoot.git/blame - TFluka/TFluka.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Correction in ion definition.
[u/mrichter/AliRoot.git] / TFluka / TFluka.cxx
CommitLineData
829fb838 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$ */
17
18//
19// Realisation of the TVirtualMC interface for the FLUKA code
20// (See official web side http://www.fluka.org/).
21//
22// This implementation makes use of the TGeo geometry modeller.
23// User configuration is via automatic generation of FLUKA input cards.
24//
25// Authors:
26// A. Fasso
27// E. Futo
28// A. Gheata
29// A. Morsch
30//
31
32#include <Riostream.h>
37b09b91 33#include <TList.h>
829fb838 34
829fb838 35#include "TFluka.h"
6d0e53bf 36#include "TFlukaIon.h"
a9ea1616 37#include "TFlukaCodes.h"
829fb838 38#include "TCallf77.h" //For the fortran calls
39#include "Fdblprc.h" //(DBLPRC) fluka common
81f1d030 40#include "Fsourcm.h" //(SOURCM) fluka common
41#include "Fgenstk.h" //(GENSTK) fluka common
829fb838 42#include "Fiounit.h" //(IOUNIT) fluka common
43#include "Fpaprop.h" //(PAPROP) fluka common
44#include "Fpart.h" //(PART) fluka common
45#include "Ftrackr.h" //(TRACKR) fluka common
46#include "Fpaprop.h" //(PAPROP) fluka common
47#include "Ffheavy.h" //(FHEAVY) fluka common
3a625972 48#include "Fopphst.h" //(OPPHST) fluka common
81f1d030 49#include "Fflkstk.h" //(FLKSTK) fluka common
07f5b33e 50#include "Fstepsz.h" //(STEPSZ) fluka common
7b203b6e 51#include "Fopphst.h" //(OPPHST) fluka common
a9ea1616 52#include "Fltclcm.h" //(LTCLCM) fluka common
f2a98602 53#include "Falldlt.h" //(ALLDLT) fluka common
829fb838 54
55#include "TVirtualMC.h"
3a625972 56#include "TMCProcess.h"
829fb838 57#include "TGeoManager.h"
58#include "TGeoMaterial.h"
59#include "TGeoMedium.h"
60#include "TFlukaMCGeometry.h"
6f5667d1 61#include "TGeoMCGeometry.h"
829fb838 62#include "TFlukaCerenkov.h"
1df5fa54 63#include "TFlukaConfigOption.h"
b496f27c 64#include "TFlukaScoringOption.h"
829fb838 65#include "TLorentzVector.h"
b496f27c 66#include "TArrayI.h"
a9ea1616 67#include "TArrayD.h"
78df7be0 68#include "TDatabasePDG.h"
4aba9d66 69#include "TStopwatch.h"
70
829fb838 71
72// Fluka methods that may be needed.
73#ifndef WIN32
74# define flukam flukam_
75# define fluka_openinp fluka_openinp_
8e5bf079 76# define fluka_openout fluka_openout_
829fb838 77# define fluka_closeinp fluka_closeinp_
78# define mcihad mcihad_
79# define mpdgha mpdgha_
2047b055 80# define newplo newplo_
4aba9d66 81# define genout genout_
82# define flkend flkend_
829fb838 83#else
84# define flukam FLUKAM
85# define fluka_openinp FLUKA_OPENINP
8e5bf079 86# define fluka_openout FLUKA_OPENOUT
829fb838 87# define fluka_closeinp FLUKA_CLOSEINP
88# define mcihad MCIHAD
89# define mpdgha MPDGHA
eea53470 90# define newplo NEWPLO
4aba9d66 91# define genout GENOUT
92# define flkend FLKEND
829fb838 93#endif
94
95extern "C"
96{
97 //
98 // Prototypes for FLUKA functions
99 //
100 void type_of_call flukam(const int&);
eea53470 101 void type_of_call newplo();
4aba9d66 102 void type_of_call genout();
103 void type_of_call flkend();
829fb838 104 void type_of_call fluka_openinp(const int&, DEFCHARA);
8e5bf079 105 void type_of_call fluka_openout(const int&, DEFCHARA);
829fb838 106 void type_of_call fluka_closeinp(const int&);
107 int type_of_call mcihad(const int&);
108 int type_of_call mpdgha(const int&);
109}
110
111//
112// Class implementation for ROOT
113//
114ClassImp(TFluka)
115
116//
117//----------------------------------------------------------------------------
118// TFluka constructors and destructors.
119//______________________________________________________________________________
120TFluka::TFluka()
121 :TVirtualMC(),
122 fVerbosityLevel(0),
4aba9d66 123 fNEvent(0),
1df5fa54 124 fInputFileName(""),
4aba9d66 125 fCoreInputFileName(""),
126 fCaller(kNoCaller),
127 fIcode(kNoProcess),
128 fNewReg(-1),
129 fRull(0),
130 fXsco(0),
131 fYsco(0),
132 fZsco(0),
01e832c7 133 fPItime(0),
134 fPIlength(0),
135 fNPI(0),
4aba9d66 136 fTrackIsEntering(kFALSE),
137 fTrackIsExiting(kFALSE),
138 fTrackIsNew(kFALSE),
139 fFieldFlag(kTRUE),
4aba9d66 140 fDummyBoundary(kFALSE),
141 fStopped(kFALSE),
142 fStopEvent(kFALSE),
143 fStopRun(kFALSE),
5125d6e5 144 fPrimaryElectronIndex(-1),
1b7bf6a6 145 fLowEnergyNeutronTransport(kFALSE),
4aba9d66 146 fMaterials(0),
147 fNVolumes(0),
148 fCurrentFlukaRegion(-1),
149 fNCerenkov(0),
150 fGeom(0),
151 fMCGeo(0),
fb2cbbec 152 fUserConfig(0),
6d0e53bf 153 fUserScore(0),
d3f271aa 154 fUserIons(0)
829fb838 155{
156 //
157 // Default constructor
158 //
82a3f706 159 for (Int_t i = 0; i < 4; i++) fPint[i] = 0.;
829fb838 160}
161
162//______________________________________________________________________________
163TFluka::TFluka(const char *title, Int_t verbosity, Bool_t isRootGeometrySupported)
164 :TVirtualMC("TFluka",title, isRootGeometrySupported),
165 fVerbosityLevel(verbosity),
4aba9d66 166 fNEvent(0),
829fb838 167 fInputFileName(""),
4aba9d66 168 fCoreInputFileName(""),
169 fCaller(kNoCaller),
170 fIcode(kNoProcess),
171 fNewReg(-1),
172 fRull(0),
173 fXsco(0),
174 fYsco(0),
175 fZsco(0),
01e832c7 176 fPItime(0),
177 fPIlength(0),
178 fNPI(0),
4aba9d66 179 fTrackIsEntering(kFALSE),
180 fTrackIsExiting(kFALSE),
181 fTrackIsNew(kFALSE),
182 fFieldFlag(kTRUE),
4aba9d66 183 fDummyBoundary(kFALSE),
184 fStopped(kFALSE),
185 fStopEvent(kFALSE),
186 fStopRun(kFALSE),
5125d6e5 187 fPrimaryElectronIndex(-1),
1b7bf6a6 188 fLowEnergyNeutronTransport(kFALSE),
4aba9d66 189 fMaterials(0),
190 fNVolumes(0),
191 fCurrentFlukaRegion(-1),
192 fNCerenkov(0),
193 fGeom(0),
194 fMCGeo(0),
fb2cbbec 195 fUserConfig(new TObjArray(100)),
6d0e53bf 196 fUserScore(new TObjArray(100)),
d3f271aa 197 fUserIons(0)
829fb838 198{
199 // create geometry interface
82a3f706 200 for (Int_t i = 0; i < 4; i++) fPint[i] = 0.;
201
7f13be31 202 if (fVerbosityLevel >=3)
203 cout << "<== TFluka::TFluka(" << title << ") constructor called." << endl;
204 SetCoreInputFileName();
205 SetInputFileName();
11e4ab84 206 fMCGeo = new TGeoMCGeometry("MCGeo", "TGeo Implementation of VirtualMCGeometry", kFALSE);
fb2cbbec 207 fGeom = new TFlukaMCGeometry("geom", "FLUKA VMC Geometry");
829fb838 208 if (verbosity > 2) fGeom->SetDebugMode(kTRUE);
8e5bf079 209 PrintHeader();
829fb838 210}
211
212//______________________________________________________________________________
4aba9d66 213TFluka::~TFluka()
214{
215 // Destructor
1df5fa54 216 if (fVerbosityLevel >=3)
4aba9d66 217 cout << "<== TFluka::~TFluka() destructor called." << endl;
218 if (fMaterials) delete [] fMaterials;
1df5fa54 219
eac7af60 220// delete fGeom;
221// delete fMCGeo;
1df5fa54 222
fb2cbbec 223 if (fUserConfig) {
4aba9d66 224 fUserConfig->Delete();
225 delete fUserConfig;
1df5fa54 226 }
6d184c54 227
228 if (fUserScore) {
4aba9d66 229 fUserScore->Delete();
230 delete fUserScore;
6d184c54 231 }
829fb838 232}
233
234//
235//______________________________________________________________________________
236// TFluka control methods
237//______________________________________________________________________________
238void TFluka::Init() {
239//
240// Geometry initialisation
241//
242 if (fVerbosityLevel >=3) cout << "==> TFluka::Init() called." << endl;
243
244 if (!gGeoManager) new TGeoManager("geom", "FLUKA geometry");
245 fApplication->ConstructGeometry();
d59acfe7 246 if (!gGeoManager->IsClosed()) {
247 TGeoVolume *top = (TGeoVolume*)gGeoManager->GetListOfVolumes()->First();
248 gGeoManager->SetTopVolume(top);
249 gGeoManager->CloseGeometry("di");
250 } else {
251 TGeoNodeCache *cache = gGeoManager->GetCache();
252 if (!cache->HasIdArray()) {
a9ea1616 253 Warning("Init", "Node ID tracking must be enabled with TFluka: enabling...\n");
d59acfe7 254 cache->BuildIdArray();
255 }
256 }
829fb838 257 fNVolumes = fGeom->NofVolumes();
258 fGeom->CreateFlukaMatFile("flukaMat.inp");
259 if (fVerbosityLevel >=3) {
260 printf("== Number of volumes: %i\n ==", fNVolumes);
261 cout << "\t* InitPhysics() - Prepare input file to be called" << endl;
6d184c54 262 }
881cb248 263
264 fApplication->InitGeometry();
661663fa 265 fApplication->ConstructOpGeometry();
78df7be0 266 //
267 // Add ions to PDG Data base
268 //
269 AddParticlesToPdgDataBase();
129dede3 270 fApplication->AddIons();
cee6a756 271 //
829fb838 272}
273
274
275//______________________________________________________________________________
276void TFluka::FinishGeometry() {
277//
278// Build-up table with region to medium correspondance
279//
280 if (fVerbosityLevel >=3) {
281 cout << "==> TFluka::FinishGeometry() called." << endl;
2753cb27 282 printf("----FinishGeometry - applying misalignment if any\n");
829fb838 283 cout << "<== TFluka::FinishGeometry() called." << endl;
284 }
2753cb27 285 TVirtualMCApplication::Instance()->MisalignGeometry();
829fb838 286}
287
288//______________________________________________________________________________
289void TFluka::BuildPhysics() {
290//
291// Prepare FLUKA input files and call FLUKA physics initialisation
292//
293
294 if (fVerbosityLevel >=3)
4aba9d66 295 cout << "==> TFluka::BuildPhysics() called." << endl;
6d184c54 296
297
298 if (fVerbosityLevel >=3) {
4aba9d66 299 TList *medlist = gGeoManager->GetListOfMedia();
300 TIter next(medlist);
301 TGeoMedium* med = 0x0;
302 TGeoMaterial* mat = 0x0;
303 Int_t ic = 0;
304
305 while((med = (TGeoMedium*)next()))
306 {
307 mat = med->GetMaterial();
308 printf("Medium %5d %12s %5d %5d\n", ic, (med->GetName()), med->GetId(), mat->GetIndex());
309 ic++;
310 }
6d184c54 311 }
312
d23f4fcd 313
6d184c54 314 // Prepare input file with the current physics settings
315
829fb838 316 InitPhysics();
b8a8a88c 317// Open fortran files
829fb838 318 const char* fname = fInputFileName;
319 fluka_openinp(lunin, PASSCHARA(fname));
8e5bf079 320 fluka_openout(11, PASSCHARA("fluka.out"));
b8a8a88c 321// Read input cards
4aba9d66 322 cout << "==> TFluka::BuildPhysics() Read input cards." << endl;
323 TStopwatch timer;
324 timer.Start();
b8a8a88c 325 GLOBAL.lfdrtr = true;
829fb838 326 flukam(1);
4aba9d66 327 cout << "<== TFluka::BuildPhysics() Read input cards End"
328 << Form(" R:%.2fs C:%.2fs", timer.RealTime(),timer.CpuTime()) << endl;
b8a8a88c 329// Close input file
829fb838 330 fluka_closeinp(lunin);
b8a8a88c 331// Finish geometry
829fb838 332 FinishGeometry();
829fb838 333}
334
335//______________________________________________________________________________
336void TFluka::ProcessEvent() {
337//
338// Process one event
339//
b496f27c 340 if (fStopRun) {
4aba9d66 341 Warning("ProcessEvent", "User Run Abortion: No more events handled !\n");
342 fNEvent += 1;
343 return;
b496f27c 344 }
345
346 if (fVerbosityLevel >=3)
4aba9d66 347 cout << "==> TFluka::ProcessEvent() called." << endl;
b496f27c 348 fApplication->GeneratePrimaries();
81f1d030 349 SOURCM.lsouit = true;
b496f27c 350 flukam(1);
351 if (fVerbosityLevel >=3)
4aba9d66 352 cout << "<== TFluka::ProcessEvent() called." << endl;
b496f27c 353 //
354 // Increase event number
355 //
356 fNEvent += 1;
829fb838 357}
358
359//______________________________________________________________________________
360Bool_t TFluka::ProcessRun(Int_t nevent) {
361//
362// Run steering
363//
4678abb9 364
829fb838 365 if (fVerbosityLevel >=3)
366 cout << "==> TFluka::ProcessRun(" << nevent << ") called."
4aba9d66 367 << endl;
829fb838 368
369 if (fVerbosityLevel >=2) {
370 cout << "\t* GLOBAL.fdrtr = " << (GLOBAL.lfdrtr?'T':'F') << endl;
371 cout << "\t* Calling flukam again..." << endl;
372 }
373
829fb838 374 Int_t todo = TMath::Abs(nevent);
375 for (Int_t ev = 0; ev < todo; ev++) {
4aba9d66 376 TStopwatch timer;
377 timer.Start();
829fb838 378 fApplication->BeginEvent();
379 ProcessEvent();
380 fApplication->FinishEvent();
4aba9d66 381 cout << "Event: "<< ev
382 << Form(" R:%.2fs C:%.2fs", timer.RealTime(),timer.CpuTime()) << endl;
829fb838 383 }
384
385 if (fVerbosityLevel >=3)
386 cout << "<== TFluka::ProcessRun(" << nevent << ") called."
4aba9d66 387 << endl;
388
eea53470 389 // Write fluka specific scoring output
4aba9d66 390 genout();
eea53470 391 newplo();
4aba9d66 392 flkend();
eea53470 393
829fb838 394 return kTRUE;
395}
396
397//_____________________________________________________________________________
398// methods for building/management of geometry
399
400// functions from GCONS
401//____________________________________________________________________________
402void TFluka::Gfmate(Int_t imat, char *name, Float_t &a, Float_t &z,
4aba9d66 403 Float_t &dens, Float_t &radl, Float_t &absl,
404 Float_t* /*ubuf*/, Int_t& /*nbuf*/) {
829fb838 405//
406 TGeoMaterial *mat;
407 TIter next (gGeoManager->GetListOfMaterials());
408 while ((mat = (TGeoMaterial*)next())) {
409 if (mat->GetUniqueID() == (UInt_t)imat) break;
410 }
411 if (!mat) {
412 Error("Gfmate", "no material with index %i found", imat);
413 return;
414 }
415 sprintf(name, "%s", mat->GetName());
416 a = mat->GetA();
417 z = mat->GetZ();
418 dens = mat->GetDensity();
419 radl = mat->GetRadLen();
420 absl = mat->GetIntLen();
421}
422
423//______________________________________________________________________________
424void TFluka::Gfmate(Int_t imat, char *name, Double_t &a, Double_t &z,
4aba9d66 425 Double_t &dens, Double_t &radl, Double_t &absl,
426 Double_t* /*ubuf*/, Int_t& /*nbuf*/) {
829fb838 427//
428 TGeoMaterial *mat;
429 TIter next (gGeoManager->GetListOfMaterials());
430 while ((mat = (TGeoMaterial*)next())) {
431 if (mat->GetUniqueID() == (UInt_t)imat) break;
432 }
433 if (!mat) {
434 Error("Gfmate", "no material with index %i found", imat);
435 return;
436 }
437 sprintf(name, "%s", mat->GetName());
438 a = mat->GetA();
439 z = mat->GetZ();
440 dens = mat->GetDensity();
441 radl = mat->GetRadLen();
442 absl = mat->GetIntLen();
443}
444
445// detector composition
446//______________________________________________________________________________
447void TFluka::Material(Int_t& kmat, const char* name, Double_t a,
4aba9d66 448 Double_t z, Double_t dens, Double_t radl, Double_t absl,
449 Float_t* buf, Int_t nwbuf) {
829fb838 450//
451 Double_t* dbuf = fGeom->CreateDoubleArray(buf, nwbuf);
452 Material(kmat, name, a, z, dens, radl, absl, dbuf, nwbuf);
453 delete [] dbuf;
454}
455
456//______________________________________________________________________________
457void TFluka::Material(Int_t& kmat, const char* name, Double_t a,
4aba9d66 458 Double_t z, Double_t dens, Double_t radl, Double_t absl,
459 Double_t* /*buf*/, Int_t /*nwbuf*/) {
829fb838 460//
fb2cbbec 461// Define a material
829fb838 462 TGeoMaterial *mat;
463 kmat = gGeoManager->GetListOfMaterials()->GetSize();
464 if ((z-Int_t(z)) > 1E-3) {
465 mat = fGeom->GetMakeWrongMaterial(z);
466 if (mat) {
467 mat->SetRadLen(radl,absl);
468 mat->SetUniqueID(kmat);
469 return;
470 }
471 }
472 gGeoManager->Material(name, a, z, dens, kmat, radl, absl);
473}
474
475//______________________________________________________________________________
476void TFluka::Mixture(Int_t& kmat, const char *name, Float_t *a,
4aba9d66 477 Float_t *z, Double_t dens, Int_t nlmat, Float_t *wmat) {
fb2cbbec 478//
479// Define a material mixture
829fb838 480//
481 Double_t* da = fGeom->CreateDoubleArray(a, TMath::Abs(nlmat));
482 Double_t* dz = fGeom->CreateDoubleArray(z, TMath::Abs(nlmat));
483 Double_t* dwmat = fGeom->CreateDoubleArray(wmat, TMath::Abs(nlmat));
484
485 Mixture(kmat, name, da, dz, dens, nlmat, dwmat);
486 for (Int_t i=0; i<nlmat; i++) {
487 a[i] = da[i]; z[i] = dz[i]; wmat[i] = dwmat[i];
488 }
489
490 delete [] da;
491 delete [] dz;
492 delete [] dwmat;
493}
494
495//______________________________________________________________________________
496void TFluka::Mixture(Int_t& kmat, const char *name, Double_t *a,
4aba9d66 497 Double_t *z, Double_t dens, Int_t nlmat, Double_t *wmat) {
829fb838 498//
499 // Defines mixture OR COMPOUND IMAT as composed by
500 // THE BASIC NLMAT materials defined by arrays A,Z and WMAT
501 //
502 // If NLMAT > 0 then wmat contains the proportion by
503 // weights of each basic material in the mixture.
504 //
505 // If nlmat < 0 then WMAT contains the number of atoms
506 // of a given kind into the molecule of the COMPOUND
507 // In this case, WMAT in output is changed to relative
508 // weigths.
509 //
510 Int_t i,j;
511 if (nlmat < 0) {
512 nlmat = - nlmat;
513 Double_t amol = 0;
514 for (i=0;i<nlmat;i++) {
515 amol += a[i]*wmat[i];
516 }
517 for (i=0;i<nlmat;i++) {
518 wmat[i] *= a[i]/amol;
519 }
520 }
521 kmat = gGeoManager->GetListOfMaterials()->GetSize();
522 // Check if we have elements with fractional Z
523 TGeoMaterial *mat = 0;
524 TGeoMixture *mix = 0;
525 Bool_t mixnew = kFALSE;
526 for (i=0; i<nlmat; i++) {
527 if (z[i]-Int_t(z[i]) < 1E-3) continue;
528 // We have found an element with fractional Z -> loop mixtures to look for it
529 for (j=0; j<kmat; j++) {
530 mat = (TGeoMaterial*)gGeoManager->GetListOfMaterials()->At(j);
531 if (!mat) break;
532 if (!mat->IsMixture()) continue;
533 mix = (TGeoMixture*)mat;
534 if (TMath::Abs(z[i]-mix->GetZ()) >1E-3) continue;
829fb838 535 mixnew = kTRUE;
536 break;
537 }
538 if (!mixnew) Warning("Mixture","%s : cannot find component %i with fractional Z=%f\n", name, i, z[i]);
539 break;
540 }
541 if (mixnew) {
542 Int_t nlmatnew = nlmat+mix->GetNelements()-1;
543 Double_t *anew = new Double_t[nlmatnew];
544 Double_t *znew = new Double_t[nlmatnew];
545 Double_t *wmatnew = new Double_t[nlmatnew];
546 Int_t ind=0;
547 for (j=0; j<nlmat; j++) {
548 if (j==i) continue;
549 anew[ind] = a[j];
550 znew[ind] = z[j];
551 wmatnew[ind] = wmat[j];
552 ind++;
553 }
554 for (j=0; j<mix->GetNelements(); j++) {
555 anew[ind] = mix->GetAmixt()[j];
556 znew[ind] = mix->GetZmixt()[j];
557 wmatnew[ind] = wmat[i]*mix->GetWmixt()[j];
558 ind++;
559 }
560 Mixture(kmat, name, anew, znew, dens, nlmatnew, wmatnew);
561 delete [] anew;
562 delete [] znew;
563 delete [] wmatnew;
564 return;
565 }
566 // Now we need to compact identical elements within the mixture
567 // First check if this happens
568 mixnew = kFALSE;
569 for (i=0; i<nlmat-1; i++) {
570 for (j=i+1; j<nlmat; j++) {
571 if (z[i] == z[j]) {
572 mixnew = kTRUE;
573 break;
574 }
575 }
576 if (mixnew) break;
577 }
578 if (mixnew) {
579 Int_t nlmatnew = 0;
580 Double_t *anew = new Double_t[nlmat];
581 Double_t *znew = new Double_t[nlmat];
582 memset(znew, 0, nlmat*sizeof(Double_t));
583 Double_t *wmatnew = new Double_t[nlmat];
584 Bool_t skipi;
585 for (i=0; i<nlmat; i++) {
586 skipi = kFALSE;
587 for (j=0; j<nlmatnew; j++) {
588 if (z[i] == z[j]) {
589 wmatnew[j] += wmat[i];
590 skipi = kTRUE;
591 break;
592 }
593 }
594 if (skipi) continue;
595 anew[nlmatnew] = a[i];
596 znew[nlmatnew] = z[i];
597 wmatnew[nlmatnew] = wmat[i];
598 nlmatnew++;
599 }
600 Mixture(kmat, name, anew, znew, dens, nlmatnew, wmatnew);
601 delete [] anew;
602 delete [] znew;
603 delete [] wmatnew;
604 return;
605 }
a8e4986c 606 gGeoManager->Mixture(name, a, z, dens, nlmat, wmat, kmat);
829fb838 607}
608
609//______________________________________________________________________________
610void TFluka::Medium(Int_t& kmed, const char *name, Int_t nmat,
4aba9d66 611 Int_t isvol, Int_t ifield, Double_t fieldm, Double_t tmaxfd,
612 Double_t stemax, Double_t deemax, Double_t epsil,
613 Double_t stmin, Float_t* ubuf, Int_t nbuf) {
b2129742 614 // Define a medium
615 //
829fb838 616 kmed = gGeoManager->GetListOfMedia()->GetSize()+1;
617 fMCGeo->Medium(kmed, name, nmat, isvol, ifield, fieldm, tmaxfd, stemax, deemax,
4aba9d66 618 epsil, stmin, ubuf, nbuf);
829fb838 619}
620
621//______________________________________________________________________________
622void TFluka::Medium(Int_t& kmed, const char *name, Int_t nmat,
4aba9d66 623 Int_t isvol, Int_t ifield, Double_t fieldm, Double_t tmaxfd,
624 Double_t stemax, Double_t deemax, Double_t epsil,
625 Double_t stmin, Double_t* ubuf, Int_t nbuf) {
b2129742 626 // Define a medium
627 //
829fb838 628 kmed = gGeoManager->GetListOfMedia()->GetSize()+1;
629 fMCGeo->Medium(kmed, name, nmat, isvol, ifield, fieldm, tmaxfd, stemax, deemax,
4aba9d66 630 epsil, stmin, ubuf, nbuf);
829fb838 631}
632
633//______________________________________________________________________________
634void TFluka::Matrix(Int_t& krot, Double_t thetaX, Double_t phiX,
4aba9d66 635 Double_t thetaY, Double_t phiY, Double_t thetaZ,
636 Double_t phiZ) {
637//
829fb838 638 krot = gGeoManager->GetListOfMatrices()->GetEntriesFast();
639 fMCGeo->Matrix(krot, thetaX, phiX, thetaY, phiY, thetaZ, phiZ);
640}
641
642//______________________________________________________________________________
643void TFluka::Gstpar(Int_t itmed, const char* param, Double_t parval) {
644//
645//
7b203b6e 646//
829fb838 647 Bool_t process = kFALSE;
acf2e119 648 Bool_t modelp = kFALSE;
649
829fb838 650 if (strncmp(param, "DCAY", 4) == 0 ||
651 strncmp(param, "PAIR", 4) == 0 ||
652 strncmp(param, "COMP", 4) == 0 ||
653 strncmp(param, "PHOT", 4) == 0 ||
654 strncmp(param, "PFIS", 4) == 0 ||
655 strncmp(param, "DRAY", 4) == 0 ||
656 strncmp(param, "ANNI", 4) == 0 ||
657 strncmp(param, "BREM", 4) == 0 ||
658 strncmp(param, "MUNU", 4) == 0 ||
659 strncmp(param, "CKOV", 4) == 0 ||
660 strncmp(param, "HADR", 4) == 0 ||
661 strncmp(param, "LOSS", 4) == 0 ||
662 strncmp(param, "MULS", 4) == 0 ||
695d8af9 663 strncmp(param, "RAYL", 4) == 0 ||
664 strncmp(param, "STRA", 4) == 0)
829fb838 665 {
666 process = kTRUE;
667 }
81f1d030 668
acf2e119 669 if (strncmp(param, "PRIMIO_N", 8) == 0 ||
670 strncmp(param, "PRIMIO_E", 8) == 0)
671 {
672 modelp = kTRUE;
673 }
674
829fb838 675 if (process) {
acf2e119 676 // Process switch
81f1d030 677 SetProcess(param, Int_t (parval), itmed);
acf2e119 678 } else if (modelp) {
679 // Model parameters
680 SetModelParameter(param, parval, itmed);
829fb838 681 } else {
acf2e119 682 // Cuts
81f1d030 683 SetCut(param, parval, itmed);
829fb838 684 }
acf2e119 685
686
829fb838 687}
688
689// functions from GGEOM
690//_____________________________________________________________________________
691void TFluka::Gsatt(const char *name, const char *att, Int_t val)
692{
6f5667d1 693 // Set visualisation attributes for one volume
829fb838 694 char vname[5];
695 fGeom->Vname(name,vname);
696 char vatt[5];
697 fGeom->Vname(att,vatt);
698 gGeoManager->SetVolumeAttribute(vname, vatt, val);
699}
700
701//______________________________________________________________________________
702Int_t TFluka::Gsvolu(const char *name, const char *shape, Int_t nmed,
4aba9d66 703 Float_t *upar, Int_t np) {
829fb838 704//
705 return fMCGeo->Gsvolu(name, shape, nmed, upar, np);
706}
707
708//______________________________________________________________________________
709Int_t TFluka::Gsvolu(const char *name, const char *shape, Int_t nmed,
4aba9d66 710 Double_t *upar, Int_t np) {
829fb838 711//
712 return fMCGeo->Gsvolu(name, shape, nmed, upar, np);
713}
714
715//______________________________________________________________________________
716void TFluka::Gsdvn(const char *name, const char *mother, Int_t ndiv,
4aba9d66 717 Int_t iaxis) {
829fb838 718//
719 fMCGeo->Gsdvn(name, mother, ndiv, iaxis);
720}
721
722//______________________________________________________________________________
723void TFluka::Gsdvn2(const char *name, const char *mother, Int_t ndiv,
4aba9d66 724 Int_t iaxis, Double_t c0i, Int_t numed) {
829fb838 725//
726 fMCGeo->Gsdvn2(name, mother, ndiv, iaxis, c0i, numed);
727}
728
729//______________________________________________________________________________
730void TFluka::Gsdvt(const char *name, const char *mother, Double_t step,
4aba9d66 731 Int_t iaxis, Int_t numed, Int_t ndvmx) {
732//
829fb838 733 fMCGeo->Gsdvt(name, mother, step, iaxis, numed, ndvmx);
734}
735
736//______________________________________________________________________________
737void TFluka::Gsdvt2(const char *name, const char *mother, Double_t step,
4aba9d66 738 Int_t iaxis, Double_t c0, Int_t numed, Int_t ndvmx) {
829fb838 739//
740 fMCGeo->Gsdvt2(name, mother, step, iaxis, c0, numed, ndvmx);
741}
742
743//______________________________________________________________________________
744void TFluka::Gsord(const char * /*name*/, Int_t /*iax*/) {
745//
746// Nothing to do with TGeo
747}
748
749//______________________________________________________________________________
750void TFluka::Gspos(const char *name, Int_t nr, const char *mother,
4aba9d66 751 Double_t x, Double_t y, Double_t z, Int_t irot,
752 const char *konly) {
829fb838 753//
754 fMCGeo->Gspos(name, nr, mother, x, y, z, irot, konly);
755}
756
757//______________________________________________________________________________
758void TFluka::Gsposp(const char *name, Int_t nr, const char *mother,
4aba9d66 759 Double_t x, Double_t y, Double_t z, Int_t irot,
760 const char *konly, Float_t *upar, Int_t np) {
829fb838 761 //
762 fMCGeo->Gsposp(name, nr, mother, x, y, z, irot, konly, upar, np);
763}
764
765//______________________________________________________________________________
766void TFluka::Gsposp(const char *name, Int_t nr, const char *mother,
4aba9d66 767 Double_t x, Double_t y, Double_t z, Int_t irot,
768 const char *konly, Double_t *upar, Int_t np) {
829fb838 769 //
770 fMCGeo->Gsposp(name, nr, mother, x, y, z, irot, konly, upar, np);
771}
772
773//______________________________________________________________________________
774void TFluka::Gsbool(const char* /*onlyVolName*/, const char* /*manyVolName*/) {
775//
776// Nothing to do with TGeo
777}
778
a9ea1616 779//______________________________________________________________________
780Bool_t TFluka::GetTransformation(const TString &volumePath,TGeoHMatrix &mat)
781{
782 // Returns the Transformation matrix between the volume specified
783 // by the path volumePath and the Top or mater volume. The format
784 // of the path volumePath is as follows (assuming ALIC is the Top volume)
785 // "/ALIC_1/DDIP_1/S05I_2/S05H_1/S05G_3". Here ALIC is the top most
786 // or master volume which has only 1 instance of. Of all of the daughter
787 // volumes of ALICE, DDIP volume copy #1 is indicated. Similarly for
788 // the daughter volume of DDIP is S05I copy #2 and so on.
789 // Inputs:
790 // TString& volumePath The volume path to the specific volume
791 // for which you want the matrix. Volume name
792 // hierarchy is separated by "/" while the
793 // copy number is appended using a "_".
794 // Outputs:
795 // TGeoHMatrix &mat A matrix with its values set to those
796 // appropriate to the Local to Master transformation
797 // Return:
798 // A logical value if kFALSE then an error occurred and no change to
799 // mat was made.
800
801 // We have to preserve the modeler state
802 return fMCGeo->GetTransformation(volumePath, mat);
803}
804
805//______________________________________________________________________
806Bool_t TFluka::GetShape(const TString &volumePath,TString &shapeType,
807 TArrayD &par)
808{
809 // Returns the shape and its parameters for the volume specified
810 // by volumeName.
811 // Inputs:
812 // TString& volumeName The volume name
813 // Outputs:
814 // TString &shapeType Shape type
815 // TArrayD &par A TArrayD of parameters with all of the
816 // parameters of the specified shape.
817 // Return:
818 // A logical indicating whether there was an error in getting this
819 // information
820 return fMCGeo->GetShape(volumePath, shapeType, par);
821}
822
823//______________________________________________________________________
824Bool_t TFluka::GetMaterial(const TString &volumeName,
825 TString &name,Int_t &imat,
826 Double_t &a,Double_t &z,Double_t &dens,
827 Double_t &radl,Double_t &inter,TArrayD &par)
828{
829 // Returns the Material and its parameters for the volume specified
830 // by volumeName.
831 // Note, Geant3 stores and uses mixtures as an element with an effective
832 // Z and A. Consequently, if the parameter Z is not integer, then
833 // this material represents some sort of mixture.
834 // Inputs:
835 // TString& volumeName The volume name
836 // Outputs:
837 // TSrting &name Material name
838 // Int_t &imat Material index number
839 // Double_t &a Average Atomic mass of material
840 // Double_t &z Average Atomic number of material
841 // Double_t &dens Density of material [g/cm^3]
842 // Double_t &radl Average radiation length of material [cm]
843 // Double_t &inter Average interaction length of material [cm]
844 // TArrayD &par A TArrayD of user defined parameters.
845 // Return:
846 // kTRUE if no errors
847 return fMCGeo->GetMaterial(volumeName,name,imat,a,z,dens,radl,inter,par);
848}
849
850//______________________________________________________________________
851Bool_t TFluka::GetMedium(const TString &volumeName,TString &name,
852 Int_t &imed,Int_t &nmat,Int_t &isvol,Int_t &ifield,
853 Double_t &fieldm,Double_t &tmaxfd,Double_t &stemax,
854 Double_t &deemax,Double_t &epsil, Double_t &stmin,
855 TArrayD &par)
856{
857 // Returns the Medium and its parameters for the volume specified
858 // by volumeName.
859 // Inputs:
860 // TString& volumeName The volume name.
861 // Outputs:
862 // TString &name Medium name
863 // Int_t &nmat Material number defined for this medium
864 // Int_t &imed The medium index number
865 // Int_t &isvol volume number defined for this medium
866 // Int_t &iflield Magnetic field flag
867 // Double_t &fieldm Magnetic field strength
868 // Double_t &tmaxfd Maximum angle of deflection per step
869 // Double_t &stemax Maximum step size
870 // Double_t &deemax Maximum fraction of energy allowed to be lost
871 // to continuous process.
872 // Double_t &epsil Boundary crossing precision
873 // Double_t &stmin Minimum step size allowed
874 // TArrayD &par A TArrayD of user parameters with all of the
875 // parameters of the specified medium.
876 // Return:
877 // kTRUE if there where no errors
878 return fMCGeo->GetMedium(volumeName,name,imed,nmat,isvol,ifield,fieldm,tmaxfd,stemax,deemax,epsil,stmin,par);
879}
880
829fb838 881//______________________________________________________________________________
882void TFluka::SetCerenkov(Int_t itmed, Int_t npckov, Float_t* ppckov,
4aba9d66 883 Float_t* absco, Float_t* effic, Float_t* rindex) {
829fb838 884//
885// Set Cerenkov properties for medium itmed
886//
887// npckov: number of sampling points
888// ppckov: energy values
889// absco: absorption length
890// effic: quantum efficiency
891// rindex: refraction index
892//
893//
894//
895// Create object holding Cerenkov properties
b6a89226 896//
897
829fb838 898 TFlukaCerenkov* cerenkovProperties = new TFlukaCerenkov(npckov, ppckov, absco, effic, rindex);
899//
900// Pass object to medium
901 TGeoMedium* medium = gGeoManager->GetMedium(itmed);
902 medium->SetCerenkovProperties(cerenkovProperties);
903}
904
b2be0e73 905void TFluka::SetCerenkov(Int_t itmed, Int_t npckov, Float_t* ppckov,
4aba9d66 906 Float_t* absco, Float_t* effic, Float_t* rindex, Float_t* rfl) {
b2be0e73 907//
908// Set Cerenkov properties for medium itmed
909//
910// npckov: number of sampling points
911// ppckov: energy values
912// absco: absorption length
913// effic: quantum efficiency
914// rindex: refraction index
915// rfl: reflectivity for boundary to medium itmed
916//
917//
918// Create object holding Cerenkov properties
919//
920 TFlukaCerenkov* cerenkovProperties = new TFlukaCerenkov(npckov, ppckov, absco, effic, rindex, rfl);
921//
922// Pass object to medium
923 TGeoMedium* medium = gGeoManager->GetMedium(itmed);
924 medium->SetCerenkovProperties(cerenkovProperties);
925}
926
927
829fb838 928//______________________________________________________________________________
b6a89226 929void TFluka::SetCerenkov(Int_t itmed, Int_t npckov, Double_t *ppckov,
930 Double_t *absco, Double_t *effic, Double_t *rindex) {
931//
932// Set Cerenkov properties for medium itmed
829fb838 933//
b6a89226 934// npckov: number of sampling points
935// ppckov: energy values
936// absco: absorption length
937// effic: quantum efficiency
938// rindex: refraction index
939//
940
941//
942// Double_t version
943 Float_t* fppckov = CreateFloatArray(ppckov, npckov);
944 Float_t* fabsco = CreateFloatArray(absco, npckov);
945 Float_t* feffic = CreateFloatArray(effic, npckov);
946 Float_t* frindex = CreateFloatArray(rindex, npckov);
947
948 SetCerenkov(itmed, npckov, fppckov, fabsco, feffic, frindex);
949
950 delete [] fppckov;
951 delete [] fabsco;
952 delete [] feffic;
953 delete [] frindex;
829fb838 954}
b2be0e73 955
b6a89226 956void TFluka::SetCerenkov(Int_t itmed, Int_t npckov, Double_t* ppckov,
957 Double_t* absco, Double_t* effic, Double_t* rindex, Double_t* rfl) {
b2be0e73 958//
b6a89226 959// Set Cerenkov properties for medium itmed
960//
961// npckov: number of sampling points
962// ppckov: energy values
963// absco: absorption length
964// effic: quantum efficiency
965// rindex: refraction index
966// rfl: reflectivity for boundary to medium itmed
967//
968
969//
970// // Double_t version
971 Float_t* fppckov = CreateFloatArray(ppckov, npckov);
972 Float_t* fabsco = CreateFloatArray(absco, npckov);
973 Float_t* feffic = CreateFloatArray(effic, npckov);
974 Float_t* frindex = CreateFloatArray(rindex, npckov);
975 Float_t* frfl = CreateFloatArray(rfl, npckov);
976
977 SetCerenkov(itmed, npckov, fppckov, fabsco, feffic, frindex, frfl);
978
979 delete [] fppckov;
980 delete [] fabsco;
981 delete [] feffic;
982 delete [] frindex;
983 delete [] frfl;
b2be0e73 984}
985
829fb838 986// Euclid
987//______________________________________________________________________________
988void TFluka::WriteEuclid(const char* /*fileName*/, const char* /*topVol*/,
989 Int_t /*number*/, Int_t /*nlevel*/) {
990//
991// Not with TGeo
a9ea1616 992 Warning("WriteEuclid", "Not implemented !");
829fb838 993}
994
995
996
997//_____________________________________________________________________________
998// methods needed by the stepping
999//____________________________________________________________________________
1000
1001Int_t TFluka::GetMedium() const {
1002//
1003// Get the medium number for the current fluka region
1004//
ab2afdda 1005 if (gGeoManager->IsOutside()) {
1006 return (-1);
1007 } else {
1008 return (fGeom->GetMedium()); // this I need to check due to remapping !!!
1009 }
829fb838 1010}
1011
a9ea1616 1012//____________________________________________________________________________
1013Int_t TFluka::GetDummyRegion() const
1014{
1015// Returns index of the dummy region.
1016 return fGeom->GetDummyRegion();
1017}
829fb838 1018
a9ea1616 1019//____________________________________________________________________________
1020Int_t TFluka::GetDummyLattice() const
1021{
1022// Returns index of the dummy lattice.
1023 return fGeom->GetDummyLattice();
1024}
829fb838 1025
1026//____________________________________________________________________________
1027// particle table usage
1028// ID <--> PDG transformations
1029//_____________________________________________________________________________
1030Int_t TFluka::IdFromPDG(Int_t pdg) const
1031{
975333ab 1032
829fb838 1033 //
1034 // Return Fluka code from PDG and pseudo ENDF code
6d0e53bf 1035 Int_t idSpecial[4] = {TFlukaIon::GetIonPdg(2,4),
1036 TFlukaIon::GetIonPdg(2,3),
1037 TFlukaIon::GetIonPdg(1,3),
1038 TFlukaIon::GetIonPdg(1,2)};
829fb838 1039 // Catch the feedback photons
a9ea1616 1040 if (pdg == 50000051) return (kFLUKAoptical);
d3f271aa 1041
1042 // Light ions (d,t,h3,alpha)
975333ab 1043 for (Int_t i = 0; i < 4; i++) {
1044 if (pdg == idSpecial[i]) return (i + kFLUKAcodemin);
1045 }
1046
d3f271aa 1047 // Heavy ions
1048 if (pdg > TFlukaIon::GetIonPdg(1,1)) return (-2);
6d0e53bf 1049
829fb838 1050 // MCIHAD() goes from pdg to fluka internal.
1051 Int_t intfluka = mcihad(pdg);
1052 // KPTOIP array goes from internal to official
1053 return GetFlukaKPTOIP(intfluka);
1054}
1055
1056//______________________________________________________________________________
1057Int_t TFluka::PDGFromId(Int_t id) const
1058{
1059 //
1060 // Return PDG code and pseudo ENDF code from Fluka code
6d0e53bf 1061 Int_t idSpecial[6] = {TFlukaIon::GetIonPdg(2,4), // alpha
1062 TFlukaIon::GetIonPdg(2,3), // He3
1063 TFlukaIon::GetIonPdg(1,3), // triton
1064 TFlukaIon::GetIonPdg(1,2), // deuteron
1065 TFlukaIon::GetIonPdg(0,0), // gen. ion
1066 50000050};
829fb838 1067 // IPTOKP array goes from official to internal
1068
a9ea1616 1069 if (id == kFLUKAoptical) {
829fb838 1070// Cerenkov photon
4aba9d66 1071// if (fVerbosityLevel >= 3)
1072// printf("\n PDGFromId: Cerenkov Photon \n");
1073 return 50000050;
829fb838 1074 }
1075// Error id
ece92b30 1076 if (id == 0 || id < kFLUKAcodemin || id > kFLUKAcodemax) {
66e5eb54 1077 if (fVerbosityLevel >= 3)
a9923346 1078 printf("PDGFromId: Error id = 0 %5d %5d\n", id, fCaller);
4aba9d66 1079 return -1;
829fb838 1080 }
1081// Good id
f926898e 1082 if (id > 0) {
4aba9d66 1083 Int_t intfluka = GetFlukaIPTOKP(id);
1084 if (intfluka == 0) {
1085 if (fVerbosityLevel >= 3)
1086 printf("PDGFromId: Error intfluka = 0: %d\n", id);
1087 return -1;
1088 } else if (intfluka < 0) {
1089 if (fVerbosityLevel >= 3)
1090 printf("PDGFromId: Error intfluka < 0: %d\n", id);
1091 return -1;
1092 }
1093// if (fVerbosityLevel >= 3)
1094// printf("mpdgha called with %d %d \n", id, intfluka);
1095 return mpdgha(intfluka);
f926898e 1096 } else {
4aba9d66 1097 // ions and optical photons
1098 return idSpecial[id - kFLUKAcodemin];
829fb838 1099 }
829fb838 1100}
1101
bd3d5c8a 1102void TFluka::StopTrack()
1103{
1104 // Set stopping conditions
1105 // Works for photons and charged particles
1106 fStopped = kTRUE;
1107}
1108
829fb838 1109//_____________________________________________________________________________
1110// methods for physics management
1111//____________________________________________________________________________
1112//
1113// set methods
1114//
1115
1df5fa54 1116void TFluka::SetProcess(const char* flagName, Int_t flagValue, Int_t imed)
829fb838 1117{
1118// Set process user flag for material imat
1df5fa54 1119//
1120//
1121// Update if already in the list
829fb838 1122//
fb2cbbec 1123 TIter next(fUserConfig);
1df5fa54 1124 TFlukaConfigOption* proc;
1125 while((proc = (TFlukaConfigOption*)next()))
1126 {
4aba9d66 1127 if (proc->Medium() == imed) {
1128 proc->SetProcess(flagName, flagValue);
1129 return;
1130 }
1df5fa54 1131 }
fb2cbbec 1132 proc = new TFlukaConfigOption(imed);
1133 proc->SetProcess(flagName, flagValue);
1134 fUserConfig->Add(proc);
1135}
1136
1137//______________________________________________________________________________
1138Bool_t TFluka::SetProcess(const char* flagName, Int_t flagValue)
1139{
1140// Set process user flag
1df5fa54 1141//
1df5fa54 1142//
fb2cbbec 1143 SetProcess(flagName, flagValue, -1);
1df5fa54 1144 return kTRUE;
829fb838 1145}
1146
1147//______________________________________________________________________________
1148void TFluka::SetCut(const char* cutName, Double_t cutValue, Int_t imed)
1149{
1150// Set user cut value for material imed
1151//
fb2cbbec 1152 TIter next(fUserConfig);
1153 TFlukaConfigOption* proc;
1154 while((proc = (TFlukaConfigOption*)next()))
1155 {
4aba9d66 1156 if (proc->Medium() == imed) {
1157 proc->SetCut(cutName, cutValue);
1158 return;
1159 }
fb2cbbec 1160 }
1161
1162 proc = new TFlukaConfigOption(imed);
1163 proc->SetCut(cutName, cutValue);
1164 fUserConfig->Add(proc);
829fb838 1165}
1166
acf2e119 1167
1168//______________________________________________________________________________
1169void TFluka::SetModelParameter(const char* parName, Double_t parValue, Int_t imed)
1170{
1171// Set model parameter for material imed
1172//
1173 TIter next(fUserConfig);
1174 TFlukaConfigOption* proc;
1175 while((proc = (TFlukaConfigOption*)next()))
1176 {
4aba9d66 1177 if (proc->Medium() == imed) {
1178 proc->SetModelParameter(parName, parValue);
1179 return;
1180 }
acf2e119 1181 }
1182
1183 proc = new TFlukaConfigOption(imed);
1184 proc->SetModelParameter(parName, parValue);
1185 fUserConfig->Add(proc);
1186}
1187
829fb838 1188//______________________________________________________________________________
1189Bool_t TFluka::SetCut(const char* cutName, Double_t cutValue)
1190{
1191// Set user cut value
1192//
1df5fa54 1193//
fb2cbbec 1194 SetCut(cutName, cutValue, -1);
1195 return kTRUE;
829fb838 1196}
1197
f450e9d0 1198
6f1aaa8e 1199void TFluka::SetUserScoring(const char* option, const char* sdum, Int_t npr, char* outfile, Float_t* what)
b496f27c 1200{
1201//
f450e9d0 1202// Adds a user scoring option to the list
b496f27c 1203//
6f1aaa8e 1204 TFlukaScoringOption* opt = new TFlukaScoringOption(option, sdum, npr,outfile,what);
f450e9d0 1205 fUserScore->Add(opt);
1206}
1207//______________________________________________________________________________
6f1aaa8e 1208void TFluka::SetUserScoring(const char* option, const char* sdum, Int_t npr, char* outfile, Float_t* what,
1209 const char* det1, const char* det2, const char* det3)
f450e9d0 1210{
1211//
1212// Adds a user scoring option to the list
1213//
6f1aaa8e 1214 TFlukaScoringOption* opt = new TFlukaScoringOption(option, sdum, npr, outfile, what, det1, det2, det3);
b496f27c 1215 fUserScore->Add(opt);
1216}
b496f27c 1217
829fb838 1218//______________________________________________________________________________
1219Double_t TFluka::Xsec(char*, Double_t, Int_t, Int_t)
1220{
a9ea1616 1221 Warning("Xsec", "Not yet implemented.!\n"); return -1.;
829fb838 1222}
1223
1224
1225//______________________________________________________________________________
1226void TFluka::InitPhysics()
1227{
1228//
1229// Physics initialisation with preparation of FLUKA input cards
1230//
fb2cbbec 1231// Construct file names
1232 FILE *pFlukaVmcCoreInp, *pFlukaVmcFlukaMat, *pFlukaVmcInp;
fb2cbbec 1233 TString sFlukaVmcTmp = "flukaMat.inp";
1234 TString sFlukaVmcInp = GetInputFileName();
ff2d1491 1235 TString sFlukaVmcCoreInp = GetCoreInputFileName();
fb2cbbec 1236
1237// Open files
1238 if ((pFlukaVmcCoreInp = fopen(sFlukaVmcCoreInp.Data(),"r")) == NULL) {
4aba9d66 1239 Warning("InitPhysics", "\nCannot open file %s\n",sFlukaVmcCoreInp.Data());
1240 exit(1);
fb2cbbec 1241 }
1242 if ((pFlukaVmcFlukaMat = fopen(sFlukaVmcTmp.Data(),"r")) == NULL) {
4aba9d66 1243 Warning("InitPhysics", "\nCannot open file %s\n",sFlukaVmcTmp.Data());
1244 exit(1);
fb2cbbec 1245 }
1246 if ((pFlukaVmcInp = fopen(sFlukaVmcInp.Data(),"w")) == NULL) {
4aba9d66 1247 Warning("InitPhysics", "\nCannot open file %s\n",sFlukaVmcInp.Data());
1248 exit(1);
fb2cbbec 1249 }
829fb838 1250
fb2cbbec 1251// Copy core input file
1252 Char_t sLine[255];
1253 Float_t fEventsPerRun;
829fb838 1254
fb2cbbec 1255 while ((fgets(sLine,255,pFlukaVmcCoreInp)) != NULL) {
4aba9d66 1256 if (strncmp(sLine,"GEOEND",6) != 0)
1257 fprintf(pFlukaVmcInp,"%s",sLine); // copy until GEOEND card
1258 else {
1259 fprintf(pFlukaVmcInp,"GEOEND\n"); // add GEOEND card
1260 goto flukamat;
1261 }
fb2cbbec 1262 } // end of while until GEOEND card
1263
829fb838 1264
fb2cbbec 1265 flukamat:
1266 while ((fgets(sLine,255,pFlukaVmcFlukaMat)) != NULL) { // copy flukaMat.inp file
4aba9d66 1267 fprintf(pFlukaVmcInp,"%s\n",sLine);
fb2cbbec 1268 }
1269
1270 while ((fgets(sLine,255,pFlukaVmcCoreInp)) != NULL) {
8fc475a1 1271 if (strncmp(sLine,"START",5) != 0)
4aba9d66 1272 fprintf(pFlukaVmcInp,"%s\n",sLine);
1273 else {
1274 sscanf(sLine+10,"%10f",&fEventsPerRun);
1275 goto fin;
1276 }
8fc475a1 1277 } //end of while until START card
fb2cbbec 1278
1279 fin:
829fb838 1280
f450e9d0 1281
1282// Pass information to configuration objects
829fb838 1283
fb2cbbec 1284 Float_t fLastMaterial = fGeom->GetLastMaterialIndex();
1285 TFlukaConfigOption::SetStaticInfo(pFlukaVmcInp, 3, fLastMaterial, fGeom);
1286
1287 TIter next(fUserConfig);
1288 TFlukaConfigOption* proc;
f450e9d0 1289 while((proc = dynamic_cast<TFlukaConfigOption*> (next()))) proc->WriteFlukaInputCards();
1290//
1291// Process Fluka specific scoring options
1292//
1293 TFlukaScoringOption::SetStaticInfo(pFlukaVmcInp, fGeom);
0bb2c369 1294 Float_t loginp = -49.0;
f450e9d0 1295 Int_t inp = 0;
1296 Int_t nscore = fUserScore->GetEntries();
1297
a9ea1616 1298 TFlukaScoringOption *mopo = 0;
1299 TFlukaScoringOption *mopi = 0;
fb2cbbec 1300
f450e9d0 1301 for (Int_t isc = 0; isc < nscore; isc++)
1302 {
4aba9d66 1303 mopo = dynamic_cast<TFlukaScoringOption*> (fUserScore->At(isc));
1304 char* fileName = mopo->GetFileName();
1305 Int_t size = strlen(fileName);
1306 Float_t lun = -1.;
f450e9d0 1307//
1308// Check if new output file has to be opened
4aba9d66 1309 for (Int_t isci = 0; isci < isc; isci++) {
1310
1311
1312 mopi = dynamic_cast<TFlukaScoringOption*> (fUserScore->At(isci));
1313 if(strncmp(mopi->GetFileName(), fileName, size)==0) {
1314 //
1315 // No, the file already exists
1316 lun = mopi->GetLun();
1317 mopo->SetLun(lun);
1318 break;
1319 }
1320 } // inner loop
1321
1322 if (lun == -1.) {
1323 // Open new output file
1324 inp++;
1325 mopo->SetLun(loginp + inp);
1326 mopo->WriteOpenFlukaFile();
1327 }
1328 mopo->WriteFlukaInputCards();
f450e9d0 1329 }
b8a8a88c 1330
1331// Add RANDOMIZ card
1332 fprintf(pFlukaVmcInp,"RANDOMIZ %10.1f%10.0f\n", 1., Float_t(gRandom->GetSeed()));
6d0e53bf 1333// User defined ion
5fc9f5d5 1334 if (fUserIons) TFlukaIon::WriteUserInputCard(pFlukaVmcInp);
8fc475a1 1335// Add START and STOP card
1336 fprintf(pFlukaVmcInp,"START %10.1f\n",fEventsPerRun);
f450e9d0 1337 fprintf(pFlukaVmcInp,"STOP \n");
829fb838 1338
1339
1340// Close files
3b8c325d 1341 fclose(pFlukaVmcCoreInp);
1342 fclose(pFlukaVmcFlukaMat);
1343 fclose(pFlukaVmcInp);
fb2cbbec 1344
1345
1346//
1347// Initialisation needed for Cerenkov photon production and transport
1348 TObjArray *matList = GetFlukaMaterials();
1349 Int_t nmaterial = matList->GetEntriesFast();
9968e86c 1350 fMaterials = new Int_t[nmaterial+25];
fb2cbbec 1351
1352 for (Int_t im = 0; im < nmaterial; im++)
1353 {
4aba9d66 1354 TGeoMaterial* material = dynamic_cast<TGeoMaterial*> (matList->At(im));
1355 Int_t idmat = material->GetIndex();
1356 fMaterials[idmat] = im;
fb2cbbec 1357 }
829fb838 1358} // end of InitPhysics
1359
1360
1361//______________________________________________________________________________
07f5b33e 1362void TFluka::SetMaxStep(Double_t step)
829fb838 1363{
07f5b33e 1364// Set the maximum step size
4aba9d66 1365// if (step > 1.e4) return;
07f5b33e 1366
4aba9d66 1367// Int_t mreg=0, latt=0;
1368// fGeom->GetCurrentRegion(mreg, latt);
ff2d1491 1369
1370
4aba9d66 1371 Int_t mreg = fGeom->GetCurrentRegion();
9c0c08ce 1372 STEPSZ.stepmx[mreg - 1] = step;
829fb838 1373}
1374
2f09b80e 1375
1376Double_t TFluka::MaxStep() const
1377{
1378// Return the maximum for current medium
1379 Int_t mreg, latt;
1380 fGeom->GetCurrentRegion(mreg, latt);
1381 return (STEPSZ.stepmx[mreg - 1]);
1382}
1383
829fb838 1384//______________________________________________________________________________
1385void TFluka::SetMaxNStep(Int_t)
1386{
1387// SetMaxNStep is dummy procedure in TFluka !
1388 if (fVerbosityLevel >=3)
1389 cout << "SetMaxNStep is dummy procedure in TFluka !" << endl;
1390}
1391
1392//______________________________________________________________________________
1393void TFluka::SetUserDecay(Int_t)
1394{
1395// SetUserDecay is dummy procedure in TFluka !
1396 if (fVerbosityLevel >=3)
1397 cout << "SetUserDecay is dummy procedure in TFluka !" << endl;
1398}
1399
1400//
1401// dynamic properties
1402//
1403//______________________________________________________________________________
1404void TFluka::TrackPosition(TLorentzVector& position) const
1405{
1406// Return the current position in the master reference frame of the
1407// track being transported
1408// TRACKR.atrack = age of the particle
1409// TRACKR.xtrack = x-position of the last point
1410// TRACKR.ytrack = y-position of the last point
1411// TRACKR.ztrack = z-position of the last point
a9ea1616 1412 FlukaCallerCode_t caller = GetCaller();
1413 if (caller == kENDRAW || caller == kUSDRAW ||
1414 caller == kBXExiting || caller == kBXEntering ||
1415 caller == kUSTCKV) {
42b936d1 1416 position.SetX(GetXsco());
1417 position.SetY(GetYsco());
1418 position.SetZ(GetZsco());
1419 position.SetT(TRACKR.atrack);
829fb838 1420 }
5125d6e5 1421 else if (caller == kMGDRAW) {
1422 Int_t i = -1;
1423 if ((i = fPrimaryElectronIndex) > -1) {
1424 // Primary Electron Ionisation
15a8a899 1425 Double_t x, y, z, t;
1426 GetPrimaryElectronPosition(i, x, y, z, t);
5125d6e5 1427 position.SetX(x);
1428 position.SetY(y);
1429 position.SetZ(z);
15a8a899 1430 position.SetT(t);
5125d6e5 1431 } else {
1432 position.SetX(TRACKR.xtrack[TRACKR.ntrack]);
1433 position.SetY(TRACKR.ytrack[TRACKR.ntrack]);
1434 position.SetZ(TRACKR.ztrack[TRACKR.ntrack]);
1435 position.SetT(TRACKR.atrack);
1436 }
829fb838 1437 }
a9ea1616 1438 else if (caller == kSODRAW) {
42b936d1 1439 Int_t ist = FLKSTK.npflka;
1440 position.SetX(FLKSTK.xflk[ist]);
1441 position.SetY(FLKSTK.yflk[ist]);
1442 position.SetZ(FLKSTK.zflk[ist]);
1443 position.SetT(FLKSTK.agestk[ist]);
a9ea1616 1444 } else if (caller == kMGResumedTrack) {
42b936d1 1445 position.SetX(TRACKR.spausr[0]);
1446 position.SetY(TRACKR.spausr[1]);
1447 position.SetZ(TRACKR.spausr[2]);
1448 position.SetT(TRACKR.spausr[3]);
829fb838 1449 }
1450 else
42b936d1 1451 Warning("TrackPosition","position not available");
829fb838 1452}
1453
1454//______________________________________________________________________________
1455void TFluka::TrackPosition(Double_t& x, Double_t& y, Double_t& z) const
1456{
1457// Return the current position in the master reference frame of the
1458// track being transported
1459// TRACKR.atrack = age of the particle
1460// TRACKR.xtrack = x-position of the last point
1461// TRACKR.ytrack = y-position of the last point
1462// TRACKR.ztrack = z-position of the last point
a9ea1616 1463 FlukaCallerCode_t caller = GetCaller();
1464 if (caller == kENDRAW || caller == kUSDRAW ||
1465 caller == kBXExiting || caller == kBXEntering ||
1466 caller == kUSTCKV) {
5125d6e5 1467 x = GetXsco();
1468 y = GetYsco();
1469 z = GetZsco();
829fb838 1470 }
42b936d1 1471 else if (caller == kMGDRAW) {
5125d6e5 1472 Int_t i = -1;
1473 if ((i = fPrimaryElectronIndex) > -1) {
15a8a899 1474 Double_t t;
1475 GetPrimaryElectronPosition(i, x, y, z, t);
5125d6e5 1476 } else {
1477 x = TRACKR.xtrack[TRACKR.ntrack];
1478 y = TRACKR.ytrack[TRACKR.ntrack];
1479 z = TRACKR.ztrack[TRACKR.ntrack];
1480 }
829fb838 1481 }
42b936d1 1482 else if (caller == kSODRAW) {
1483 Int_t ist = FLKSTK.npflka;
1484 x = FLKSTK.xflk[ist];
1485 y = FLKSTK.yflk[ist];
1486 z = FLKSTK.zflk[ist];
1487 }
a9ea1616 1488 else if (caller == kMGResumedTrack) {
42b936d1 1489 x = TRACKR.spausr[0];
1490 y = TRACKR.spausr[1];
1491 z = TRACKR.spausr[2];
5d80a015 1492 }
829fb838 1493 else
42b936d1 1494 Warning("TrackPosition","position not available");
829fb838 1495}
1496
1497//______________________________________________________________________________
1498void TFluka::TrackMomentum(TLorentzVector& momentum) const
1499{
1500// Return the direction and the momentum (GeV/c) of the track
1501// currently being transported
1502// TRACKR.ptrack = momentum of the particle (not always defined, if
1503// < 0 must be obtained from etrack)
1504// TRACKR.cx,y,ztrck = direction cosines of the current particle
1505// TRACKR.etrack = total energy of the particle
1506// TRACKR.jtrack = identity number of the particle
1507// PAPROP.am[TRACKR.jtrack] = particle mass in gev
a9ea1616 1508 FlukaCallerCode_t caller = GetCaller();
1509 FlukaProcessCode_t icode = GetIcode();
1510
82a3f706 1511 if (caller != kEEDRAW &&
1512 caller != kMGResumedTrack &&
1513 caller != kSODRAW &&
1514 caller != kUSDRAW &&
a9ea1616 1515 (caller != kENDRAW || (icode != kEMFSCOstopping1 && icode != kEMFSCOstopping2))) {
42b936d1 1516 if (TRACKR.ptrack >= 0) {
1517 momentum.SetPx(TRACKR.ptrack*TRACKR.cxtrck);
1518 momentum.SetPy(TRACKR.ptrack*TRACKR.cytrck);
1519 momentum.SetPz(TRACKR.ptrack*TRACKR.cztrck);
1520 momentum.SetE(TRACKR.etrack);
1521 return;
1522 }
1523 else {
1524 Double_t p = sqrt(TRACKR.etrack * TRACKR.etrack - ParticleMassFPC(TRACKR.jtrack) * ParticleMassFPC(TRACKR.jtrack));
1525 momentum.SetPx(p*TRACKR.cxtrck);
1526 momentum.SetPy(p*TRACKR.cytrck);
1527 momentum.SetPz(p*TRACKR.cztrck);
1528 momentum.SetE(TRACKR.etrack);
1529 return;
1530 }
a9ea1616 1531 } else if (caller == kMGResumedTrack) {
42b936d1 1532 momentum.SetPx(TRACKR.spausr[4]);
1533 momentum.SetPy(TRACKR.spausr[5]);
1534 momentum.SetPz(TRACKR.spausr[6]);
1535 momentum.SetE (TRACKR.spausr[7]);
1536 return;
a9ea1616 1537 } else if (caller == kENDRAW && (icode == kEMFSCOstopping1 || icode == kEMFSCOstopping2)) {
1538 momentum.SetPx(0.);
1539 momentum.SetPy(0.);
1540 momentum.SetPz(0.);
1541 momentum.SetE(TrackMass());
42b936d1 1542
1543 } else if (caller == kSODRAW) {
1544 Int_t ist = FLKSTK.npflka;
1545 Double_t p = FLKSTK.pmoflk[ist];
1546 Int_t ifl = FLKSTK.iloflk[ist];
1547 Double_t m = PAPROP.am[ifl + 6];
1548 Double_t e = TMath::Sqrt(p * p + m * m);
1549 momentum.SetPx(p * FLKSTK.txflk[ist]);
1550 momentum.SetPy(p * FLKSTK.tyflk[ist]);
1551 momentum.SetPz(p * FLKSTK.tzflk[ist]);
1552 momentum.SetE(e);
82a3f706 1553 } else if (caller == kUSDRAW) {
6df2c1cc 1554 if (icode == kEMFSCObrems ||
1555 icode == kEMFSCOmoller ||
1556 icode == kEMFSCObhabha ||
1557 icode == kEMFSCOcompton )
1558 {
82a3f706 1559 momentum.SetPx(fPint[0]);
1560 momentum.SetPy(fPint[1]);
1561 momentum.SetPz(fPint[2]);
1562 momentum.SetE(fPint[3]);
6df2c1cc 1563 } else if (icode == kKASKADdray ||
1564 icode == kKASKADbrems ||
1565 icode == kKASKADpair) {
1566 momentum.SetPx(GENSTK.plr[0] * GENSTK.cxr[0]);
1567 momentum.SetPy(GENSTK.plr[0] * GENSTK.cyr[0]);
1568 momentum.SetPz(GENSTK.plr[0] * GENSTK.czr[0]);
1569 momentum.SetE (GENSTK.tki[0] + PAPROP.am[GENSTK.kpart[0]+6]);
82a3f706 1570 } else {
6df2c1cc 1571 Double_t p = sqrt(TRACKR.etrack * TRACKR.etrack
1572 - ParticleMassFPC(TRACKR.jtrack)
1573 * ParticleMassFPC(TRACKR.jtrack));
82a3f706 1574 momentum.SetPx(p*TRACKR.cxtrck);
1575 momentum.SetPy(p*TRACKR.cytrck);
1576 momentum.SetPz(p*TRACKR.cztrck);
1577 momentum.SetE(TRACKR.etrack);
1578 }
829fb838 1579 }
1580 else
1581 Warning("TrackMomentum","momentum not available");
1582}
1583
1584//______________________________________________________________________________
1585void TFluka::TrackMomentum(Double_t& px, Double_t& py, Double_t& pz, Double_t& e) const
1586{
1587// Return the direction and the momentum (GeV/c) of the track
1588// currently being transported
1589// TRACKR.ptrack = momentum of the particle (not always defined, if
1590// < 0 must be obtained from etrack)
1591// TRACKR.cx,y,ztrck = direction cosines of the current particle
1592// TRACKR.etrack = total energy of the particle
1593// TRACKR.jtrack = identity number of the particle
1594// PAPROP.am[TRACKR.jtrack] = particle mass in gev
a9ea1616 1595 FlukaCallerCode_t caller = GetCaller();
1596 FlukaProcessCode_t icode = GetIcode();
42b936d1 1597 if (caller != kEEDRAW &&
1598 caller != kMGResumedTrack &&
1599 caller != kSODRAW &&
82a3f706 1600 caller != kUSDRAW &&
a9ea1616 1601 (caller != kENDRAW || (icode != kEMFSCOstopping1 && icode != kEMFSCOstopping2))) {
829fb838 1602 if (TRACKR.ptrack >= 0) {
1603 px = TRACKR.ptrack*TRACKR.cxtrck;
1604 py = TRACKR.ptrack*TRACKR.cytrck;
1605 pz = TRACKR.ptrack*TRACKR.cztrck;
a9ea1616 1606 e = TRACKR.etrack;
829fb838 1607 return;
1608 }
1609 else {
ece92b30 1610 Double_t p = sqrt(TRACKR.etrack * TRACKR.etrack - ParticleMassFPC(TRACKR.jtrack) * ParticleMassFPC(TRACKR.jtrack));
829fb838 1611 px = p*TRACKR.cxtrck;
1612 py = p*TRACKR.cytrck;
1613 pz = p*TRACKR.cztrck;
a9ea1616 1614 e = TRACKR.etrack;
829fb838 1615 return;
1616 }
a9ea1616 1617 } else if (caller == kMGResumedTrack) {
5d80a015 1618 px = TRACKR.spausr[4];
1619 py = TRACKR.spausr[5];
1620 pz = TRACKR.spausr[6];
1621 e = TRACKR.spausr[7];
0773d0ac 1622 return;
a9ea1616 1623 } else if (caller == kENDRAW && (icode == kEMFSCOstopping1 || icode == kEMFSCOstopping2)) {
1624 px = 0.;
1625 py = 0.;
1626 pz = 0.;
1627 e = TrackMass();
42b936d1 1628 } else if (caller == kSODRAW) {
1629 Int_t ist = FLKSTK.npflka;
1630 Double_t p = FLKSTK.pmoflk[ist];
1631 Int_t ifl = FLKSTK.iloflk[ist];
1632 Double_t m = PAPROP.am[ifl + 6];
1633 e = TMath::Sqrt(p * p + m * m);
1634 px = p * FLKSTK.txflk[ist];
1635 py = p * FLKSTK.tyflk[ist];
1636 pz = p * FLKSTK.tzflk[ist];
82a3f706 1637 } else if (caller == kUSDRAW) {
6df2c1cc 1638 if (icode == kEMFSCObrems ||
1639 icode == kEMFSCOmoller ||
1640 icode == kEMFSCObhabha ||
1641 icode == kEMFSCOcompton )
1642 {
82a3f706 1643 px = fPint[0];
1644 py = fPint[1];
1645 pz = fPint[2];
1646 e = fPint[3];
6df2c1cc 1647 } else if (icode == kKASKADdray ||
1648 icode == kKASKADbrems ||
1649 icode == kKASKADpair) {
1650 px = GENSTK.plr[0] * GENSTK.cxr[0];
1651 py = GENSTK.plr[0] * GENSTK.cyr[0];
1652 pz = GENSTK.plr[0] * GENSTK.czr[0];
1653 e = GENSTK.tki[0] + PAPROP.am[GENSTK.kpart[0]+6];
82a3f706 1654 } else {
1655 Double_t p = sqrt(TRACKR.etrack * TRACKR.etrack - ParticleMassFPC(TRACKR.jtrack) * ParticleMassFPC(TRACKR.jtrack));
1656 px = p*TRACKR.cxtrck;
1657 py = p*TRACKR.cytrck;
1658 pz = p*TRACKR.cztrck;
1659 e = TRACKR.etrack;
1660 }
829fb838 1661 }
1662 else
42b936d1 1663 Warning("TrackMomentum","momentum not available");
829fb838 1664}
1665
1666//______________________________________________________________________________
1667Double_t TFluka::TrackStep() const
1668{
1669// Return the length in centimeters of the current step
1670// TRACKR.ctrack = total curved path
42b936d1 1671 FlukaCallerCode_t caller = GetCaller();
01e832c7 1672 if (caller == kMGDRAW) {
1673 Int_t i;
1674 if ((i = fPrimaryElectronIndex) > -1) {
1675 if (i > 0) {
1676 return (fPIlength[i] - fPIlength[i-1]);
1677 } else {
1678 Double_t s (TRACKR.ctrack - (fPIlength[fNPI - 1] - fPIlength[0]));
1679 return s;
1680 }
1681 } else {
1682 return TRACKR.ctrack;
1683 }
1684 } else if (caller == kBXEntering || caller == kBXExiting ||
1685 caller == kENDRAW || caller == kUSDRAW ||
1686 caller == kUSTCKV || caller == kMGResumedTrack ||
1687 caller == kSODRAW)
1688 {
42b936d1 1689 return 0.0;
01e832c7 1690 } else {
1691 Warning("TrackStep", "track step not available");
1692 return 0.0;
1693 }
829fb838 1694}
1695
1696//______________________________________________________________________________
1697Double_t TFluka::TrackLength() const
1698{
1699// TRACKR.cmtrck = cumulative curved path since particle birth
a9ea1616 1700 FlukaCallerCode_t caller = GetCaller();
01e832c7 1701 if (caller == kMGDRAW) {
1702 Int_t i;
1703 if ((i = fPrimaryElectronIndex) > -1) {
1704 return fPIlength[i];
1705 } else {
1706 return TRACKR.cmtrck;
1707 }
1708
1709 } else if (caller == kBXEntering || caller == kBXExiting ||
1710 caller == kENDRAW || caller == kUSDRAW || caller == kUSTCKV)
1711 return TRACKR.cmtrck;
a9ea1616 1712 else if (caller == kMGResumedTrack)
01e832c7 1713 return TRACKR.spausr[8];
82a3f706 1714 else if (caller == kSODRAW)
1715 return 0.0;
669cede4 1716 else {
01e832c7 1717 Warning("TrackLength", "track length not available for caller %5d \n", caller);
1718 return 0.0;
1719 }
829fb838 1720}
1721
01e832c7 1722
829fb838 1723//______________________________________________________________________________
1724Double_t TFluka::TrackTime() const
1725{
1726// Return the current time of flight of the track being transported
1727// TRACKR.atrack = age of the particle
a9ea1616 1728 FlukaCallerCode_t caller = GetCaller();
15a8a899 1729 if (caller == kMGDRAW) {
1730 Int_t i;
1731 if ((i = fPrimaryElectronIndex) > -1) {
01e832c7 1732 Double_t t = fPItime[i];
15a8a899 1733 return t;
1734 } else {
1735 return TRACKR.atrack;
1736 }
1737 } else if (caller == kBXEntering || caller == kBXExiting ||
1738 caller == kENDRAW || caller == kUSDRAW ||
1739 caller == kUSTCKV)
829fb838 1740 return TRACKR.atrack;
a9ea1616 1741 else if (caller == kMGResumedTrack)
5d80a015 1742 return TRACKR.spausr[3];
42b936d1 1743 else if (caller == kSODRAW) {
1744 return (FLKSTK.agestk[FLKSTK.npflka]);
1745 }
669cede4 1746 else {
1747 Warning("TrackTime", "track time not available");
1748 return 0.0;
1749 }
829fb838 1750}
1751
1752//______________________________________________________________________________
1753Double_t TFluka::Edep() const
1754{
1755// Energy deposition
1756// if TRACKR.ntrack = 0, TRACKR.mtrack = 0:
1757// -->local energy deposition (the value and the point are not recorded in TRACKR)
1758// but in the variable "rull" of the procedure "endraw.cxx"
1759// if TRACKR.ntrack > 0, TRACKR.mtrack = 0:
1760// -->no energy loss along the track
1761// if TRACKR.ntrack > 0, TRACKR.mtrack > 0:
1762// -->energy loss distributed along the track
07f5b33e 1763// TRACKR.dtrack = energy deposition of the jth deposition event
829fb838 1764
1765 // If coming from bxdraw we have 2 steps of 0 length and 0 edep
669cede4 1766 // If coming from usdraw we just signal particle production - no edep
1767 // If just first time after resuming, no edep for the primary
a9ea1616 1768 FlukaCallerCode_t caller = GetCaller();
ada781c7 1769
a9ea1616 1770 if (caller == kBXExiting || caller == kBXEntering ||
42b936d1 1771 caller == kUSDRAW || caller == kMGResumedTrack ||
1772 caller == kSODRAW)
1773 return 0.0;
829fb838 1774 Double_t sum = 0;
5125d6e5 1775 Int_t i = -1;
09cdde8a 1776
ada781c7 1777 // Material with primary ionisation activated but number of primary electrons nprim = 0
1778 if (fPrimaryElectronIndex == -2) return 0.0;
1779 // nprim > 0
5125d6e5 1780 if ((i = fPrimaryElectronIndex) > -1) {
1781 // Primary ionisation
ada781c7 1782 sum = GetPrimaryElectronKineticEnergy(i);
1783 if (sum > 100.) {
1784 printf("edep > 100. %d %d %f \n", i, ALLDLT.nalldl, sum);
1785 }
1786 return sum;
5125d6e5 1787 } else {
1788 // Normal ionisation
1789 if (TRACKR.mtrack > 1) printf("Edep: %6d\n", TRACKR.mtrack);
1790
1791 for ( Int_t j=0;j<TRACKR.mtrack;j++) {
1792 sum +=TRACKR.dtrack[j];
1793 }
1794 if (TRACKR.ntrack == 0 && TRACKR.mtrack == 0)
1795 return fRull + sum;
1796 else {
1797 return sum;
1798 }
829fb838 1799 }
1800}
1801
18e0cabb 1802//______________________________________________________________________________
1803Int_t TFluka::CorrectFlukaId() const
1804{
1805 // since we don't put photons and e- created bellow transport cut on the vmc stack
1806 // and there is a call to endraw for energy deposition for each of them
1807 // and they have the track number of their parent, but different identity (pdg)
4aba9d66 1808 // so we want to assign also their parent identity.
cc7af78a 1809
a9923346 1810 if( (IsTrackStop())
18e0cabb 1811 && TRACKR.ispusr[mkbmx2 - 4] == TRACKR.ispusr[mkbmx2 - 1]
1812 && TRACKR.jtrack != TRACKR.ispusr[mkbmx2 - 3] ) {
1813 if (fVerbosityLevel >=3)
1814 cout << "CorrectFlukaId() for icode=" << GetIcode()
1815 << " track=" << TRACKR.ispusr[mkbmx2 - 1]
1816 << " current PDG=" << PDGFromId(TRACKR.jtrack)
1817 << " assign parent PDG=" << PDGFromId(TRACKR.ispusr[mkbmx2 - 3]) << endl;
1818 return TRACKR.ispusr[mkbmx2 - 3]; // assign parent identity
1819 }
13858fbd 1820 if (TRACKR.jtrack <= 64){
cc7af78a 1821 return TRACKR.jtrack;
1822 } else {
1823 return TRACKR.j0trck;
1824 }
18e0cabb 1825}
1826
1827
829fb838 1828//______________________________________________________________________________
1829Int_t TFluka::TrackPid() const
1830{
1831// Return the id of the particle transported
1832// TRACKR.jtrack = identity number of the particle
a9ea1616 1833 FlukaCallerCode_t caller = GetCaller();
42b936d1 1834 if (caller != kEEDRAW && caller != kSODRAW) {
18e0cabb 1835 return PDGFromId( CorrectFlukaId() );
f926898e 1836 }
42b936d1 1837 else if (caller == kSODRAW) {
1838 return PDGFromId(FLKSTK.iloflk[FLKSTK.npflka]);
1839 }
829fb838 1840 else
1841 return -1000;
1842}
1843
1844//______________________________________________________________________________
1845Double_t TFluka::TrackCharge() const
1846{
1847// Return charge of the track currently transported
1848// PAPROP.ichrge = electric charge of the particle
1849// TRACKR.jtrack = identity number of the particle
13858fbd 1850
a9ea1616 1851 FlukaCallerCode_t caller = GetCaller();
42b936d1 1852 if (caller != kEEDRAW && caller != kSODRAW)
1853 return PAPROP.ichrge[CorrectFlukaId() + 6];
1854 else if (caller == kSODRAW) {
1855 Int_t ifl = PDGFromId(FLKSTK.iloflk[FLKSTK.npflka]);
1856 return PAPROP.ichrge[ifl + 6];
1857 }
829fb838 1858 else
1859 return -1000.0;
1860}
1861
1862//______________________________________________________________________________
1863Double_t TFluka::TrackMass() const
1864{
1865// PAPROP.am = particle mass in GeV
1866// TRACKR.jtrack = identity number of the particle
a9ea1616 1867 FlukaCallerCode_t caller = GetCaller();
42b936d1 1868 if (caller != kEEDRAW && caller != kSODRAW)
18e0cabb 1869 return PAPROP.am[CorrectFlukaId()+6];
42b936d1 1870 else if (caller == kSODRAW) {
82a3f706 1871 Int_t ifl = FLKSTK.iloflk[FLKSTK.npflka];
42b936d1 1872 return PAPROP.am[ifl + 6];
1873 }
829fb838 1874 else
1875 return -1000.0;
1876}
1877
1878//______________________________________________________________________________
1879Double_t TFluka::Etot() const
1880{
1881// TRACKR.etrack = total energy of the particle
6df2c1cc 1882 FlukaCallerCode_t caller = GetCaller();
1883 FlukaProcessCode_t icode = GetIcode();
1884 if (caller != kEEDRAW && caller != kSODRAW && caller != kUSDRAW)
1885 {
1886 return TRACKR.etrack;
1887 } else if (caller == kUSDRAW) {
1888 if (icode == kEMFSCObrems ||
1889 icode == kEMFSCOmoller ||
1890 icode == kEMFSCObhabha ||
1891 icode == kEMFSCOcompton ) {
1892 return fPint[3];
1893 }
1894 else if (icode == kKASKADdray ||
1895 icode == kKASKADbrems ||
1896 icode == kKASKADpair) {
1897 return (GENSTK.tki[0] + PAPROP.am[GENSTK.kpart[0]+6]);
ca01d0af 1898 } else {
1899 return TRACKR.etrack;
6df2c1cc 1900 }
ca01d0af 1901
6df2c1cc 1902 }
42b936d1 1903 else if (caller == kSODRAW) {
1904 Int_t ist = FLKSTK.npflka;
1905 Double_t p = FLKSTK.pmoflk[ist];
1906 Int_t ifl = FLKSTK.iloflk[ist];
1907 Double_t m = PAPROP.am[ifl + 6];
1908 Double_t e = TMath::Sqrt(p * p + m * m);
1909 return e;
1910 }
ca01d0af 1911 printf("Etot %5d %5d \n", caller, icode);
6df2c1cc 1912
1913 return -1000.0;
829fb838 1914}
1915
1916//
1917// track status
1918//
1919//______________________________________________________________________________
1920Bool_t TFluka::IsNewTrack() const
1921{
1922// Return true for the first call of Stepping()
1923 return fTrackIsNew;
1924}
1925
0dabe425 1926void TFluka::SetTrackIsNew(Bool_t flag)
1927{
1928// Return true for the first call of Stepping()
1929 fTrackIsNew = flag;
1930
1931}
1932
1933
829fb838 1934//______________________________________________________________________________
1935Bool_t TFluka::IsTrackInside() const
1936{
1937// True if the track is not at the boundary of the current volume
1938// In Fluka a step is always inside one kind of material
1939// If the step would go behind the region of one material,
1940// it will be shortened to reach only the boundary.
1941// Therefore IsTrackInside() is always true.
a9ea1616 1942 FlukaCallerCode_t caller = GetCaller();
1943 if (caller == kBXEntering || caller == kBXExiting)
829fb838 1944 return 0;
1945 else
1946 return 1;
1947}
1948
1949//______________________________________________________________________________
1950Bool_t TFluka::IsTrackEntering() const
1951{
1952// True if this is the first step of the track in the current volume
1953
a9ea1616 1954 FlukaCallerCode_t caller = GetCaller();
1955 if (caller == kBXEntering)
829fb838 1956 return 1;
1957 else return 0;
1958}
1959
1960//______________________________________________________________________________
1961Bool_t TFluka::IsTrackExiting() const
1962{
1963// True if track is exiting volume
1964//
a9ea1616 1965 FlukaCallerCode_t caller = GetCaller();
1966 if (caller == kBXExiting)
829fb838 1967 return 1;
1968 else return 0;
1969}
1970
1971//______________________________________________________________________________
1972Bool_t TFluka::IsTrackOut() const
1973{
1974// True if the track is out of the setup
1975// means escape
a9ea1616 1976 FlukaProcessCode_t icode = GetIcode();
1977
1978 if (icode == kKASKADescape ||
1979 icode == kEMFSCOescape ||
1980 icode == kKASNEUescape ||
1981 icode == kKASHEAescape ||
1982 icode == kKASOPHescape)
1983 return 1;
829fb838 1984 else return 0;
1985}
1986
1987//______________________________________________________________________________
1988Bool_t TFluka::IsTrackDisappeared() const
1989{
a9ea1616 1990// All inelastic interactions and decays
829fb838 1991// fIcode from usdraw
a9ea1616 1992 FlukaProcessCode_t icode = GetIcode();
1993 if (icode == kKASKADinelint || // inelastic interaction
1994 icode == kKASKADdecay || // particle decay
1995 icode == kKASKADdray || // delta ray generation by hadron
1996 icode == kKASKADpair || // direct pair production
1997 icode == kKASKADbrems || // bremsstrahlung (muon)
1998 icode == kEMFSCObrems || // bremsstrahlung (electron)
1999 icode == kEMFSCOmoller || // Moller scattering
2000 icode == kEMFSCObhabha || // Bhaba scattering
2001 icode == kEMFSCOanniflight || // in-flight annihilation
2002 icode == kEMFSCOannirest || // annihilation at rest
2003 icode == kEMFSCOpair || // pair production
2004 icode == kEMFSCOcompton || // Compton scattering
2005 icode == kEMFSCOphotoel || // Photoelectric effect
2006 icode == kKASNEUhadronic || // hadronic interaction
2047b055 2007 icode == kKASHEAdray // delta-ray
0dabe425 2008 ) return 1;
829fb838 2009 else return 0;
2010}
2011
2012//______________________________________________________________________________
2013Bool_t TFluka::IsTrackStop() const
2014{
2015// True if the track energy has fallen below the threshold
2016// means stopped by signal or below energy threshold
a9ea1616 2017 FlukaProcessCode_t icode = GetIcode();
18e0cabb 2018 if (icode == kKASKADstopping || // stopping particle
2019 icode == kKASKADtimekill || // time kill
2020 icode == kEMFSCOstopping1 || // below user-defined cut-off
2021 icode == kEMFSCOstopping2 || // below user cut-off
2022 icode == kEMFSCOtimekill || // time kill
2023 icode == kKASNEUstopping || // neutron below threshold
2024 icode == kKASNEUtimekill || // time kill
2025 icode == kKASHEAtimekill || // time kill
2026 icode == kKASOPHtimekill) return 1; // time kill
829fb838 2027 else return 0;
2028}
2029
2030//______________________________________________________________________________
2031Bool_t TFluka::IsTrackAlive() const
2032{
695d8af9 2033// Means not disappeared or not out
2034 FlukaProcessCode_t icode = GetIcode();
2035
2036 if (IsTrackOut() ||
2037 IsTrackStop() ||
2038 icode == kKASKADinelint || // inelastic interaction
2039 icode == kKASKADdecay || // particle decay
2040 icode == kEMFSCOanniflight || // in-flight annihilation
2041 icode == kEMFSCOannirest || // annihilation at rest
2042 icode == kEMFSCOpair || // pair production
2043 icode == kEMFSCOphotoel || // Photoelectric effect
2044 icode == kKASNEUhadronic // hadronic interaction
2045 )
2046 {
2047 // Exclude the cases for which the particle has disappeared (paused) but will reappear later (= alive).
2048 return 0;
2049 } else {
2050 return 1;
2051 }
829fb838 2052}
2053
2054//
2055// secondaries
2056//
2057
2058//______________________________________________________________________________
2059Int_t TFluka::NSecondaries() const
2060
2061{
2062// Number of secondary particles generated in the current step
81f1d030 2063// GENSTK.np = number of secondaries except light and heavy ions
829fb838 2064// FHEAVY.npheav = number of secondaries for light and heavy secondary ions
a9ea1616 2065 FlukaCallerCode_t caller = GetCaller();
2066 if (caller == kUSDRAW) // valid only after usdraw
4aba9d66 2067 return GENSTK.np + FHEAVY.npheav;
a9ea1616 2068 else if (caller == kUSTCKV) {
4aba9d66 2069 // Cerenkov Photon production
2070 return fNCerenkov;
7b203b6e 2071 }
829fb838 2072 return 0;
2073} // end of NSecondaries
2074
2075//______________________________________________________________________________
2076void TFluka::GetSecondary(Int_t isec, Int_t& particleId,
4aba9d66 2077 TLorentzVector& position, TLorentzVector& momentum)
829fb838 2078{
2079// Copy particles from secondary stack to vmc stack
2080//
2081
a9ea1616 2082 FlukaCallerCode_t caller = GetCaller();
2083 if (caller == kUSDRAW) { // valid only after usdraw
4aba9d66 2084 if (GENSTK.np > 0) {
2085 // Hadronic interaction
2086 if (isec >= 0 && isec < GENSTK.np) {
2087 particleId = PDGFromId(GENSTK.kpart[isec]);
2088 position.SetX(fXsco);
2089 position.SetY(fYsco);
2090 position.SetZ(fZsco);
2091 position.SetT(TRACKR.atrack);
2092 momentum.SetPx(GENSTK.plr[isec]*GENSTK.cxr[isec]);
2093 momentum.SetPy(GENSTK.plr[isec]*GENSTK.cyr[isec]);
2094 momentum.SetPz(GENSTK.plr[isec]*GENSTK.czr[isec]);
2095 momentum.SetE(GENSTK.tki[isec] + PAPROP.am[GENSTK.kpart[isec]+6]);
2096 }
2097 else if (isec >= GENSTK.np && isec < GENSTK.np + FHEAVY.npheav) {
2098 Int_t jsec = isec - GENSTK.np;
2099 particleId = FHEAVY.kheavy[jsec]; // this is Fluka id !!!
2100 position.SetX(fXsco);
2101 position.SetY(fYsco);
2102 position.SetZ(fZsco);
2103 position.SetT(TRACKR.atrack);
2104 momentum.SetPx(FHEAVY.pheavy[jsec]*FHEAVY.cxheav[jsec]);
2105 momentum.SetPy(FHEAVY.pheavy[jsec]*FHEAVY.cyheav[jsec]);
2106 momentum.SetPz(FHEAVY.pheavy[jsec]*FHEAVY.czheav[jsec]);
2107 if (FHEAVY.tkheav[jsec] >= 3 && FHEAVY.tkheav[jsec] <= 6)
2108 momentum.SetE(FHEAVY.tkheav[jsec] + PAPROP.am[jsec+6]);
2109 else if (FHEAVY.tkheav[jsec] > 6)
2110 momentum.SetE(FHEAVY.tkheav[jsec] + FHEAVY.amnhea[jsec]); // to be checked !!!
2111 }
2112 else
2113 Warning("GetSecondary","isec out of range");
2114 }
a9ea1616 2115 } else if (caller == kUSTCKV) {
4aba9d66 2116 Int_t index = OPPHST.lstopp - isec;
2117 position.SetX(OPPHST.xoptph[index]);
2118 position.SetY(OPPHST.yoptph[index]);
2119 position.SetZ(OPPHST.zoptph[index]);
2120 position.SetT(OPPHST.agopph[index]);
2121 Double_t p = OPPHST.poptph[index];
2122
2123 momentum.SetPx(p * OPPHST.txopph[index]);
2124 momentum.SetPy(p * OPPHST.tyopph[index]);
2125 momentum.SetPz(p * OPPHST.tzopph[index]);
2126 momentum.SetE(p);
829fb838 2127 }
2128 else
4aba9d66 2129 Warning("GetSecondary","no secondaries available");
7b203b6e 2130
829fb838 2131} // end of GetSecondary
2132
7b203b6e 2133
829fb838 2134//______________________________________________________________________________
2135TMCProcess TFluka::ProdProcess(Int_t) const
2136
2137{
2138// Name of the process that has produced the secondary particles
2139// in the current step
0dabe425 2140
a9ea1616 2141 Int_t mugamma = (TRACKR.jtrack == kFLUKAphoton ||
4aba9d66 2142 TRACKR.jtrack == kFLUKAmuplus ||
2143 TRACKR.jtrack == kFLUKAmuminus);
a9ea1616 2144 FlukaProcessCode_t icode = GetIcode();
2145
2146 if (icode == kKASKADdecay) return kPDecay;
2147 else if (icode == kKASKADpair || icode == kEMFSCOpair) return kPPair;
2148 else if (icode == kEMFSCOcompton) return kPCompton;
2149 else if (icode == kEMFSCOphotoel) return kPPhotoelectric;
2150 else if (icode == kKASKADbrems || icode == kEMFSCObrems) return kPBrem;
2151 else if (icode == kKASKADdray || icode == kKASHEAdray) return kPDeltaRay;
2152 else if (icode == kEMFSCOmoller || icode == kEMFSCObhabha) return kPDeltaRay;
2153 else if (icode == kEMFSCOanniflight || icode == kEMFSCOannirest) return kPAnnihilation;
2154 else if (icode == kKASKADinelint) {
4aba9d66 2155 if (!mugamma) return kPHadronic;
2156 else if (TRACKR.jtrack == kFLUKAphoton) return kPPhotoFission;
2157 else return kPMuonNuclear;
829fb838 2158 }
a9ea1616 2159 else if (icode == kEMFSCOrayleigh) return kPRayleigh;
829fb838 2160// Fluka codes 100, 300 and 400 still to be investigasted
a9ea1616 2161 else return kPNoProcess;
829fb838 2162}
2163
829fb838 2164
b496f27c 2165Int_t TFluka::StepProcesses(TArrayI &proc) const
2166{
2167 //
2168 // Return processes active in the current step
2169 //
e71bcde8 2170 FlukaProcessCode_t icode = GetIcode();
2171 FlukaCallerCode_t caller = GetCaller();
b496f27c 2172 proc.Set(1);
2173 TMCProcess iproc;
87ad3c3e 2174 if (caller == kBXEntering || caller == kBXExiting || caller == kEEDRAW || caller == kSODRAW) {
e71bcde8 2175 iproc = kPTransportation;
87ad3c3e 2176 }
2177 else if (caller == kUSTCKV) {
2178 iproc = kPCerenkov;
e71bcde8 2179 } else {
2180 switch (icode) {
2181 case kEMFSCO:
82a3f706 2182 if (Edep() > 0.) {
2183 iproc = kPEnergyLoss;
2184 } else {
2185 iproc = kPTransportation;
2186 }
e71bcde8 2187 break;
6df2c1cc 2188 case kKASKAD:
2189 if (Edep() > 0.) {
2190 iproc = kPEnergyLoss;
2191 } else {
2192 iproc = kPTransportation;
2193 }
2194 break;
87ad3c3e 2195 case kKASHEA:
2196 case kKASNEU:
2197 case kKASOPH:
2198 case kKASKADescape:
2199 case kEMFSCOescape:
2200 case kKASNEUescape:
2201 case kKASHEAescape:
2202 case kKASOPHescape:
2203 iproc = kPTransportation;
2204 break;
e71bcde8 2205 case kKASKADtimekill:
2206 case kEMFSCOtimekill:
2207 case kKASNEUtimekill:
2208 case kKASHEAtimekill:
2209 case kKASOPHtimekill:
2210 iproc = kPTOFlimit;
2211 break;
2212 case kKASKADstopping:
e71bcde8 2213 case kEMFSCOstopping1:
2214 case kEMFSCOstopping2:
e71bcde8 2215 case kKASNEUstopping:
e71bcde8 2216 iproc = kPStop;
ca01d0af 2217 break;
2218 case kKASKADinelint:
2219 case kKASNEUhadronic:
2220 iproc = kPHadronic;
2221 break;
2222 case kKASKADinelarecoil:
2223 iproc = kPHadronic;
2224 break;
2225 case kKASKADnelint:
2226 iproc = kPHElastic;
e71bcde8 2227 break;
2228 case kKASOPHabsorption:
2229 iproc = kPLightAbsorption;
2230 break;
2231 case kKASOPHrefraction:
2232 iproc = kPLightRefraction;
87ad3c3e 2233 break;
e71bcde8 2234 case kEMFSCOlocaldep :
2235 iproc = kPPhotoelectric;
2236 break;
2237 default:
2238 iproc = ProdProcess(0);
2239 }
b496f27c 2240 }
e71bcde8 2241
07f5b33e 2242 proc[0] = iproc;
b496f27c 2243 return 1;
2244}
829fb838 2245//______________________________________________________________________________
2246Int_t TFluka::VolId2Mate(Int_t id) const
2247{
2248//
2249// Returns the material number for a given volume ID
2250//
2251 return fMCGeo->VolId2Mate(id);
2252}
2253
2254//______________________________________________________________________________
2255const char* TFluka::VolName(Int_t id) const
2256{
2257//
2258// Returns the volume name for a given volume ID
2259//
2260 return fMCGeo->VolName(id);
2261}
2262
a8e4986c 2263Int_t TFluka::MediumId(const Text_t* mediumName) const
2264{
2265 //
2266 // Return the unique medium id for medium with name mediumName
2267 TList *medlist = gGeoManager->GetListOfMedia();
2268 TGeoMedium* med = (TGeoMedium*) medlist->FindObject(mediumName);
2269 if (med) {
2270 return (med->GetId());
2271 } else {
2272 return (-1);
2273 }
2274}
2275
829fb838 2276//______________________________________________________________________________
2277Int_t TFluka::VolId(const Text_t* volName) const
2278{
2279//
2280// Converts from volume name to volume ID.
2281// Time consuming. (Only used during set-up)
2282// Could be replaced by hash-table
2283//
09cd6497 2284 char sname[20];
2285 Int_t len;
2286 strncpy(sname, volName, len = strlen(volName));
2287 sname[len] = 0;
2288 while (sname[len - 1] == ' ') sname[--len] = 0;
2289 return fMCGeo->VolId(sname);
829fb838 2290}
2291
2292//______________________________________________________________________________
2293Int_t TFluka::CurrentVolID(Int_t& copyNo) const
2294{
2295//
2296// Return the logical id and copy number corresponding to the current fluka region
2297//
2298 if (gGeoManager->IsOutside()) return 0;
2299 TGeoNode *node = gGeoManager->GetCurrentNode();
2300 copyNo = node->GetNumber();
2301 Int_t id = node->GetVolume()->GetNumber();
2302 return id;
2303}
2304
2305//______________________________________________________________________________
2306Int_t TFluka::CurrentVolOffID(Int_t off, Int_t& copyNo) const
2307{
2308//
2309// Return the logical id and copy number of off'th mother
2310// corresponding to the current fluka region
2311//
2312 if (off<0 || off>gGeoManager->GetLevel()) return 0;
2313 if (off==0) return CurrentVolID(copyNo);
2314 TGeoNode *node = gGeoManager->GetMother(off);
2315 if (!node) return 0;
2316 copyNo = node->GetNumber();
2317 return node->GetVolume()->GetNumber();
2318}
2319
2320//______________________________________________________________________________
2321const char* TFluka::CurrentVolName() const
2322{
2323//
2324// Return the current volume name
2325//
ab2afdda 2326 if (gGeoManager->IsOutside()) return "OutOfWorld";
829fb838 2327 return gGeoManager->GetCurrentVolume()->GetName();
2328}
2329
2330//______________________________________________________________________________
2331const char* TFluka::CurrentVolOffName(Int_t off) const
2332{
2333//
2334// Return the volume name of the off'th mother of the current volume
2335//
2336 if (off<0 || off>gGeoManager->GetLevel()) return 0;
2337 if (off==0) return CurrentVolName();
2338 TGeoNode *node = gGeoManager->GetMother(off);
2339 if (!node) return 0;
2340 return node->GetVolume()->GetName();
2341}
2342
d59acfe7 2343const char* TFluka::CurrentVolPath() {
2344 // Return the current volume path
2345 return gGeoManager->GetPath();
2346}
829fb838 2347//______________________________________________________________________________
a60813de 2348Int_t TFluka::CurrentMaterial(Float_t & a, Float_t & z,
4aba9d66 2349 Float_t & dens, Float_t & radl, Float_t & absl) const
829fb838 2350{
2351//
a60813de 2352// Return the current medium number and material properties
829fb838 2353//
2354 Int_t copy;
2355 Int_t id = TFluka::CurrentVolID(copy);
2356 Int_t med = TFluka::VolId2Mate(id);
a60813de 2357 TGeoVolume* vol = gGeoManager->GetCurrentVolume();
2358 TGeoMaterial* mat = vol->GetMaterial();
2359 a = mat->GetA();
2360 z = mat->GetZ();
2361 dens = mat->GetDensity();
2362 radl = mat->GetRadLen();
2363 absl = mat->GetIntLen();
2364
829fb838 2365 return med;
2366}
2367
2368//______________________________________________________________________________
2369void TFluka::Gmtod(Float_t* xm, Float_t* xd, Int_t iflag)
2370{
2371// Transforms a position from the world reference frame
2372// to the current volume reference frame.
2373//
2374// Geant3 desription:
2375// ==================
2376// Computes coordinates XD (in DRS)
2377// from known coordinates XM in MRS
2378// The local reference system can be initialized by
2379// - the tracking routines and GMTOD used in GUSTEP
2380// - a call to GMEDIA(XM,NUMED)
2381// - a call to GLVOLU(NLEVEL,NAMES,NUMBER,IER)
2382// (inverse routine is GDTOM)
2383//
2384// If IFLAG=1 convert coordinates
2385// IFLAG=2 convert direction cosinus
2386//
2387// ---
2388 Double_t xmL[3], xdL[3];
2389 Int_t i;
2390 for (i=0;i<3;i++) xmL[i]=xm[i];
2391 if (iflag == 1) gGeoManager->MasterToLocal(xmL,xdL);
2392 else gGeoManager->MasterToLocalVect(xmL,xdL);
2393 for (i=0;i<3;i++) xd[i] = xdL[i];
2394}
2395
2396//______________________________________________________________________________
2397void TFluka::Gmtod(Double_t* xm, Double_t* xd, Int_t iflag)
2398{
2047b055 2399//
2400// See Gmtod(Float_t*, Float_t*, Int_t)
2401//
829fb838 2402 if (iflag == 1) gGeoManager->MasterToLocal(xm,xd);
2403 else gGeoManager->MasterToLocalVect(xm,xd);
2404}
2405
2406//______________________________________________________________________________
2407void TFluka::Gdtom(Float_t* xd, Float_t* xm, Int_t iflag)
2408{
2409// Transforms a position from the current volume reference frame
2410// to the world reference frame.
2411//
2412// Geant3 desription:
2413// ==================
2414// Computes coordinates XM (Master Reference System
2415// knowing the coordinates XD (Detector Ref System)
2416// The local reference system can be initialized by
2417// - the tracking routines and GDTOM used in GUSTEP
2418// - a call to GSCMED(NLEVEL,NAMES,NUMBER)
2419// (inverse routine is GMTOD)
2420//
2421// If IFLAG=1 convert coordinates
2422// IFLAG=2 convert direction cosinus
2423//
2424// ---
2425 Double_t xmL[3], xdL[3];
2426 Int_t i;
2427 for (i=0;i<3;i++) xdL[i] = xd[i];
2428 if (iflag == 1) gGeoManager->LocalToMaster(xdL,xmL);
2429 else gGeoManager->LocalToMasterVect(xdL,xmL);
2430 for (i=0;i<3;i++) xm[i]=xmL[i];
2431}
2432
2433//______________________________________________________________________________
2434void TFluka::Gdtom(Double_t* xd, Double_t* xm, Int_t iflag)
2435{
2047b055 2436//
2437// See Gdtom(Float_t*, Float_t*, Int_t)
2438//
829fb838 2439 if (iflag == 1) gGeoManager->LocalToMaster(xd,xm);
2440 else gGeoManager->LocalToMasterVect(xd,xm);
2441}
2442
2443//______________________________________________________________________________
2444TObjArray *TFluka::GetFlukaMaterials()
2445{
2047b055 2446//
2447// Get array of Fluka materials
829fb838 2448 return fGeom->GetMatList();
2449}
2450
2451//______________________________________________________________________________
a9ea1616 2452void TFluka::SetMreg(Int_t l, Int_t lttc)
829fb838 2453{
2454// Set current fluka region
2455 fCurrentFlukaRegion = l;
a9ea1616 2456 fGeom->SetMreg(l,lttc);
829fb838 2457}
2458
2459
b496f27c 2460
2461
4aba9d66 2462//______________________________________________________________________________
b496f27c 2463TString TFluka::ParticleName(Int_t pdg) const
2464{
2465 // Return particle name for particle with pdg code pdg.
2466 Int_t ifluka = IdFromPDG(pdg);
ece92b30 2467 return TString((CHPPRP.btype[ifluka - kFLUKAcodemin]), 8);
b496f27c 2468}
2469
2470
4aba9d66 2471//______________________________________________________________________________
b496f27c 2472Double_t TFluka::ParticleMass(Int_t pdg) const
2473{
2474 // Return particle mass for particle with pdg code pdg.
2475 Int_t ifluka = IdFromPDG(pdg);
ece92b30 2476 return (PAPROP.am[ifluka - kFLUKAcodemin]);
2477}
2478
4aba9d66 2479//______________________________________________________________________________
ece92b30 2480Double_t TFluka::ParticleMassFPC(Int_t fpc) const
2481{
2482 // Return particle mass for particle with Fluka particle code fpc
2483 return (PAPROP.am[fpc - kFLUKAcodemin]);
b496f27c 2484}
2485
4aba9d66 2486//______________________________________________________________________________
b496f27c 2487Double_t TFluka::ParticleCharge(Int_t pdg) const
2488{
2489 // Return particle charge for particle with pdg code pdg.
2490 Int_t ifluka = IdFromPDG(pdg);
ece92b30 2491 return Double_t(PAPROP.ichrge[ifluka - kFLUKAcodemin]);
b496f27c 2492}
2493
4aba9d66 2494//______________________________________________________________________________
b496f27c 2495Double_t TFluka::ParticleLifeTime(Int_t pdg) const
2496{
2497 // Return particle lifetime for particle with pdg code pdg.
2498 Int_t ifluka = IdFromPDG(pdg);
ece92b30 2499 return (PAPROP.tmnlf[ifluka - kFLUKAcodemin]);
b496f27c 2500}
2501
4aba9d66 2502//______________________________________________________________________________
b496f27c 2503void TFluka::Gfpart(Int_t pdg, char* name, Int_t& type, Float_t& mass, Float_t& charge, Float_t& tlife)
2504{
2505 // Retrieve particle properties for particle with pdg code pdg.
2506
2507 strcpy(name, ParticleName(pdg).Data());
2508 type = ParticleMCType(pdg);
2509 mass = ParticleMass(pdg);
2510 charge = ParticleCharge(pdg);
2511 tlife = ParticleLifeTime(pdg);
2512}
2513
4aba9d66 2514//______________________________________________________________________________
8e5bf079 2515void TFluka::PrintHeader()
2516{
2517 //
2518 // Print a header
2519 printf("\n");
2520 printf("\n");
2521 printf("------------------------------------------------------------------------------\n");
2522 printf("- You are using the TFluka Virtual Monte Carlo Interface to FLUKA. -\n");
2523 printf("- Please see the file fluka.out for FLUKA output and licensing information. -\n");
2524 printf("------------------------------------------------------------------------------\n");
2525 printf("\n");
2526 printf("\n");
2527}
2528
b496f27c 2529
81f1d030 2530#define pshckp pshckp_
2531#define ustckv ustckv_
3a625972 2532
2533
2534extern "C" {
81f1d030 2535 void pshckp(Double_t & px, Double_t & py, Double_t & pz, Double_t & e,
4aba9d66 2536 Double_t & vx, Double_t & vy, Double_t & vz, Double_t & tof,
2537 Double_t & polx, Double_t & poly, Double_t & polz, Double_t & wgt, Int_t& ntr)
81f1d030 2538 {
2539 //
2540 // Pushes one cerenkov photon to the stack
2541 //
2542
2543 TFluka* fluka = (TFluka*) gMC;
2544 TVirtualMCStack* cppstack = fluka->GetStack();
2545 Int_t parent = TRACKR.ispusr[mkbmx2-1];
2546 cppstack->PushTrack(0, parent, 50000050,
4aba9d66 2547 px, py, pz, e,
2548 vx, vy, vz, tof,
2549 polx, poly, polz,
2550 kPCerenkov, ntr, wgt, 0);
2551 if (fluka->GetVerbosityLevel() >= 3)
2552 printf("pshckp: track=%d parent=%d lattc=%d %s\n", ntr, parent, TRACKR.lt1trk, fluka->CurrentVolName());
81f1d030 2553 }
2554
2555 void ustckv(Int_t & nphot, Int_t & mreg, Double_t & x, Double_t & y, Double_t & z)
7b203b6e 2556 {
4aba9d66 2557 //
2558 // Calls stepping in order to signal cerenkov production
2559 //
2560 TFluka *fluka = (TFluka*)gMC;
2561 fluka->SetMreg(mreg, TRACKR.lt1trk); //LTCLCM.mlatm1);
2562 fluka->SetXsco(x);
2563 fluka->SetYsco(y);
2564 fluka->SetZsco(z);
2565 fluka->SetNCerenkov(nphot);
2566 fluka->SetCaller(kUSTCKV);
2567 if (fluka->GetVerbosityLevel() >= 3)
2568 printf("ustckv: %10d mreg=%d lattc=%d newlat=%d (%f, %f, %f) edep=%f vol=%s\n",
2569 nphot, mreg, TRACKR.lt1trk, LTCLCM.newlat, x, y, z, fluka->Edep(), fluka->CurrentVolName());
2570
2571 // check region lattice consistency (debug Ernesto)
2572 // *****************************************************
2573 Int_t nodeId;
2574 Int_t volId = fluka->CurrentVolID(nodeId);
2575 Int_t crtlttc = gGeoManager->GetCurrentNodeId()+1;
2576
2577 if( mreg != volId && !gGeoManager->IsOutside() ) {
2578 cout << " ustckv: track=" << TRACKR.ispusr[mkbmx2-1] << " pdg=" << fluka->PDGFromId(TRACKR.jtrack)
2579 << " icode=" << fluka->GetIcode() << " gNstep=" << fluka->GetNstep() << endl
2580 << " fluka mreg=" << mreg << " mlttc=" << TRACKR.lt1trk << endl
2581 << " TGeo volId=" << volId << " crtlttc=" << crtlttc << endl
2582 << " common TRACKR lt1trk=" << TRACKR.lt1trk << " lt2trk=" << TRACKR.lt2trk << endl
2583 << " common LTCLCM newlat=" << LTCLCM.newlat << " mlatld=" << LTCLCM.mlatld << endl
2584 << " mlatm1=" << LTCLCM.mlatm1 << " mltsen=" << LTCLCM.mltsen << endl
2585 << " mltsm1=" << LTCLCM.mltsm1 << " mlattc=" << LTCLCM.mlattc << endl;
2586 if( TRACKR.lt1trk == crtlttc ) cout << " *************************************************************" << endl;
2587 }
2588 // *****************************************************
2589
2590
2591
2592 (TVirtualMCApplication::Instance())->Stepping();
7b203b6e 2593 }
3a625972 2594}
a9ea1616 2595
4aba9d66 2596//______________________________________________________________________________
78df7be0 2597void TFluka::AddParticlesToPdgDataBase() const
2598{
2599
2600//
2601// Add particles to the PDG data base
2602
2603 TDatabasePDG *pdgDB = TDatabasePDG::Instance();
2604
78df7be0 2605 const Double_t kAu2Gev = 0.9314943228;
2606 const Double_t khSlash = 1.0545726663e-27;
2607 const Double_t kErg2Gev = 1/1.6021773349e-3;
2608 const Double_t khShGev = khSlash*kErg2Gev;
2609 const Double_t kYear2Sec = 3600*24*365.25;
2610//
2611// Ions
2612//
78df7be0 2613 pdgDB->AddParticle("Deuteron","Deuteron",2*kAu2Gev+8.071e-3,kTRUE,
6d0e53bf 2614 0,3,"Ion",TFlukaIon::GetIonPdg(1,2));
78df7be0 2615 pdgDB->AddParticle("Triton","Triton",3*kAu2Gev+14.931e-3,kFALSE,
6d0e53bf 2616 khShGev/(12.33*kYear2Sec),3,"Ion",TFlukaIon::GetIonPdg(1,3));
78df7be0 2617 pdgDB->AddParticle("Alpha","Alpha",4*kAu2Gev+2.424e-3,kTRUE,
6d0e53bf 2618 khShGev/(12.33*kYear2Sec),6,"Ion",TFlukaIon::GetIonPdg(2,4));
78df7be0 2619 pdgDB->AddParticle("HE3","HE3",3*kAu2Gev+14.931e-3,kFALSE,
6d0e53bf 2620 0,6,"Ion",TFlukaIon::GetIonPdg(2,3));
6d0e53bf 2621
cee6a756 2622//
2623//
2624//
2625// Special particles
2626//
2627 pdgDB->AddParticle("Cherenkov","Cherenkov",0,kFALSE,
2628 0,0,"Special",GetSpecialPdg(50));
2629 pdgDB->AddParticle("FeedbackPhoton","FeedbackPhoton",0,kFALSE,
2630 0,0,"Special",GetSpecialPdg(51));
78df7be0 2631}
2632
ca01d0af 2633
4aba9d66 2634//
2635// Info about primary ionization electrons
2636//
2637
2638//______________________________________________________________________________
2639Int_t TFluka::GetNPrimaryElectrons()
f2a98602 2640{
2641 // Get number of primary electrons
2642 return ALLDLT.nalldl;
2643}
2644
4aba9d66 2645//______________________________________________________________________________
5125d6e5 2646Double_t TFluka::GetPrimaryElectronKineticEnergy(Int_t i) const
f2a98602 2647{
2648 // Returns kinetic energy of primary electron i
5125d6e5 2649
2650 Double_t ekin = -1.;
ea262cc6 2651
f2a98602 2652 if (i >= 0 && i < ALLDLT.nalldl) {
6c854012 2653 ekin = ALLDLT.talldl[i];
f2a98602 2654 } else {
4aba9d66 2655 Warning("GetPrimaryElectronKineticEnergy",
2656 "Primary electron index out of range %d %d \n",
2657 i, ALLDLT.nalldl);
f2a98602 2658 }
f0734960 2659 return ekin;
f2a98602 2660}
5125d6e5 2661
15a8a899 2662void TFluka::GetPrimaryElectronPosition(Int_t i, Double_t& x, Double_t& y, Double_t& z, Double_t& t) const
5125d6e5 2663{
2664 // Returns position of primary electron i
2665 if (i >= 0 && i < ALLDLT.nalldl) {
6c854012 2666 x = ALLDLT.xalldl[i];
2667 y = ALLDLT.yalldl[i];
2668 z = ALLDLT.zalldl[i];
15a8a899 2669 t = ALLDLT.talldl[i];
5125d6e5 2670 return;
2671 } else {
2672 Warning("GetPrimaryElectronPosition",
2673 "Primary electron index out of range %d %d \n",
2674 i, ALLDLT.nalldl);
2675 return;
2676 }
2677 return;
2678}
2679
cee6a756 2680
2681//__________________________________________________________________
2682Int_t TFluka::GetSpecialPdg(Int_t number) const
2683{
2684// Numbering for special particles
2685
2686 return 50000000 + number;
2687}
2688
13858fbd 2689
ea262cc6 2690void TFluka::PrimaryIonisationStepping(Int_t nprim)
2691{
2692// Call Stepping for primary ionisation electrons
ea262cc6 2693// Protection against nprim > mxalld
ea262cc6 2694// Multiple steps for nprim > 0
01e832c7 2695 Int_t i;
2696 fNPI = nprim;
ea262cc6 2697 if (nprim > 0) {
01e832c7 2698 CalcPrimaryIonisationTime();
ea262cc6 2699 for (i = 0; i < nprim; i++) {
2700 SetCurrentPrimaryElectronIndex(i);
2701 (TVirtualMCApplication::Instance())->Stepping();
2702 if (i == 0) SetTrackIsNew(kFALSE);
2703 }
2704 } else {
2705 // No primary electron ionisation
2706 // Call Stepping anyway but flag nprim = 0 as index = -2
2707 SetCurrentPrimaryElectronIndex(-2);
2708 (TVirtualMCApplication::Instance())->Stepping();
2709 }
2710 // Reset the index
2711 SetCurrentPrimaryElectronIndex(-1);
2712}
b6a89226 2713
2714//______________________________________________________________________
2715Float_t* TFluka::CreateFloatArray(Double_t* array, Int_t size) const
2716{
2717// Converts Double_t* array to Float_t*,
2718// !! The new array has to be deleted by user.
2719// ---
01e832c7 2720
b6a89226 2721 Float_t* floatArray;
2722 if (size>0) {
2723 floatArray = new Float_t[size];
2724 for (Int_t i=0; i<size; i++)
2725 if (array[i] >= FLT_MAX )
2726 floatArray[i] = FLT_MAX/100.;
2727 else
2728 floatArray[i] = array[i];
2729 }
2730 else {
2731 //floatArray = 0;
2732 floatArray = new Float_t[1];
2733 }
2734 return floatArray;
2735}
01e832c7 2736
2737void TFluka::CalcPrimaryIonisationTime()
2738{
2739 // Calculates the primary ionisation time
2740 if (fPItime) delete [] fPItime;
2741 fPItime = new Double_t[fNPI];
2742 if (fPIlength) delete [] fPIlength;
2743 fPIlength = new Double_t[fNPI];
2744 //
2745 Double_t px, py, pz, e, t;
2746 TrackMomentum(px, py, pz, e);
2747 Double_t p = TMath::Sqrt(px * px + py * py + pz * pz);
2748 Double_t beta = p / e;
2749 Double_t x0, y0, z0;
2750 fPItime[fNPI -1] = TRACKR.atrack;
2751 fPIlength[fNPI -1] = TRACKR.cmtrck;
2752 GetPrimaryElectronPosition(fNPI - 1, x0, y0, z0, t);
2753 if (fNPI > 1) {
2754 for (Int_t i = fNPI - 2; i > -1; i--) {
2755 Double_t x, y, z, t;
2756 GetPrimaryElectronPosition(i, x, y, z, t);
2757 Double_t ds = TMath::Sqrt((x-x0) * (x-x0) + (y-y0) * (y-y0) + (z-z0) * (z-z0));
2758 fPItime[i] = fPItime[i+1] - ds / (beta * 2.99792458e10);
2759 fPIlength[i] = fPIlength[i+1] - ds;
2760 x0 = x; y0 = y; z0 = z;
2761 }
2762 }
2763
2764}
6d0e53bf 2765
2766Bool_t TFluka::DefineIon(const char* name , Int_t z, Int_t a, Int_t q, Double_t exE, Double_t mass)
2767{
2768 // User defined ion that can be used as a primary
d3f271aa 2769 fUserIons = kTRUE;
2770 TFlukaIon::AddIon(name, z, a, q,exE, mass);
2771 return kTRUE;
6d0e53bf 2772}