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