]> 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
- Two debug routines removed.
[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>
33
829fb838 34#include "TFluka.h"
a9ea1616 35#include "TFlukaCodes.h"
829fb838 36#include "TCallf77.h" //For the fortran calls
37#include "Fdblprc.h" //(DBLPRC) fluka common
81f1d030 38#include "Fsourcm.h" //(SOURCM) fluka common
39#include "Fgenstk.h" //(GENSTK) fluka common
829fb838 40#include "Fiounit.h" //(IOUNIT) fluka common
41#include "Fpaprop.h" //(PAPROP) fluka common
42#include "Fpart.h" //(PART) fluka common
43#include "Ftrackr.h" //(TRACKR) fluka common
44#include "Fpaprop.h" //(PAPROP) fluka common
45#include "Ffheavy.h" //(FHEAVY) fluka common
3a625972 46#include "Fopphst.h" //(OPPHST) fluka common
81f1d030 47#include "Fflkstk.h" //(FLKSTK) fluka common
07f5b33e 48#include "Fstepsz.h" //(STEPSZ) fluka common
7b203b6e 49#include "Fopphst.h" //(OPPHST) fluka common
a9ea1616 50#include "Fltclcm.h" //(LTCLCM) fluka common
829fb838 51
52#include "TVirtualMC.h"
3a625972 53#include "TMCProcess.h"
829fb838 54#include "TGeoManager.h"
55#include "TGeoMaterial.h"
56#include "TGeoMedium.h"
57#include "TFlukaMCGeometry.h"
6f5667d1 58#include "TGeoMCGeometry.h"
829fb838 59#include "TFlukaCerenkov.h"
1df5fa54 60#include "TFlukaConfigOption.h"
b496f27c 61#include "TFlukaScoringOption.h"
829fb838 62#include "TLorentzVector.h"
b496f27c 63#include "TArrayI.h"
a9ea1616 64#include "TArrayD.h"
829fb838 65
66// Fluka methods that may be needed.
67#ifndef WIN32
68# define flukam flukam_
69# define fluka_openinp fluka_openinp_
8e5bf079 70# define fluka_openout fluka_openout_
829fb838 71# define fluka_closeinp fluka_closeinp_
72# define mcihad mcihad_
73# define mpdgha mpdgha_
74#else
75# define flukam FLUKAM
76# define fluka_openinp FLUKA_OPENINP
8e5bf079 77# define fluka_openout FLUKA_OPENOUT
829fb838 78# define fluka_closeinp FLUKA_CLOSEINP
79# define mcihad MCIHAD
80# define mpdgha MPDGHA
eea53470 81# define newplo NEWPLO
829fb838 82#endif
83
84extern "C"
85{
86 //
87 // Prototypes for FLUKA functions
88 //
89 void type_of_call flukam(const int&);
eea53470 90 void type_of_call newplo();
829fb838 91 void type_of_call fluka_openinp(const int&, DEFCHARA);
8e5bf079 92 void type_of_call fluka_openout(const int&, DEFCHARA);
829fb838 93 void type_of_call fluka_closeinp(const int&);
94 int type_of_call mcihad(const int&);
95 int type_of_call mpdgha(const int&);
96}
97
98//
99// Class implementation for ROOT
100//
101ClassImp(TFluka)
102
103//
104//----------------------------------------------------------------------------
105// TFluka constructors and destructors.
106//______________________________________________________________________________
107TFluka::TFluka()
108 :TVirtualMC(),
109 fVerbosityLevel(0),
1df5fa54 110 fInputFileName(""),
fb2cbbec 111 fUserConfig(0),
1df5fa54 112 fUserScore(0)
829fb838 113{
114 //
115 // Default constructor
116 //
117 fGeneratePemf = kFALSE;
118 fNVolumes = 0;
119 fCurrentFlukaRegion = -1;
a9ea1616 120 fNewReg = -1;
829fb838 121 fGeom = 0;
122 fMCGeo = 0;
123 fMaterials = 0;
124 fDummyBoundary = 0;
125 fFieldFlag = 1;
bd3d5c8a 126 fStopped = 0;
b496f27c 127 fStopEvent = 0;
128 fStopRun = 0;
129 fNEvent = 0;
829fb838 130}
131
132//______________________________________________________________________________
133TFluka::TFluka(const char *title, Int_t verbosity, Bool_t isRootGeometrySupported)
134 :TVirtualMC("TFluka",title, isRootGeometrySupported),
135 fVerbosityLevel(verbosity),
136 fInputFileName(""),
137 fTrackIsEntering(0),
138 fTrackIsExiting(0),
1df5fa54 139 fTrackIsNew(0),
fb2cbbec 140 fUserConfig(new TObjArray(100)),
1df5fa54 141 fUserScore(new TObjArray(100))
829fb838 142{
143 // create geometry interface
7f13be31 144 if (fVerbosityLevel >=3)
145 cout << "<== TFluka::TFluka(" << title << ") constructor called." << endl;
146 SetCoreInputFileName();
147 SetInputFileName();
148 SetGeneratePemf(kFALSE);
829fb838 149 fNVolumes = 0;
150 fCurrentFlukaRegion = -1;
a9ea1616 151 fNewReg = -1;
829fb838 152 fDummyBoundary = 0;
153 fFieldFlag = 1;
154 fGeneratePemf = kFALSE;
155 fMCGeo = new TGeoMCGeometry("MCGeo", "TGeo Implementation of VirtualMCGeometry", kTRUE);
fb2cbbec 156 fGeom = new TFlukaMCGeometry("geom", "FLUKA VMC Geometry");
829fb838 157 if (verbosity > 2) fGeom->SetDebugMode(kTRUE);
158 fMaterials = 0;
bd3d5c8a 159 fStopped = 0;
b496f27c 160 fStopEvent = 0;
161 fStopRun = 0;
162 fNEvent = 0;
8e5bf079 163 PrintHeader();
829fb838 164}
165
166//______________________________________________________________________________
167TFluka::~TFluka() {
168// Destructor
1df5fa54 169 if (fVerbosityLevel >=3)
170 cout << "<== TFluka::~TFluka() destructor called." << endl;
171
172 delete fGeom;
173 delete fMCGeo;
174
fb2cbbec 175 if (fUserConfig) {
176 fUserConfig->Delete();
177 delete fUserConfig;
1df5fa54 178 }
6d184c54 179
180 if (fUserScore) {
181 fUserScore->Delete();
182 delete fUserScore;
183 }
829fb838 184}
185
186//
187//______________________________________________________________________________
188// TFluka control methods
189//______________________________________________________________________________
190void TFluka::Init() {
191//
192// Geometry initialisation
193//
194 if (fVerbosityLevel >=3) cout << "==> TFluka::Init() called." << endl;
195
196 if (!gGeoManager) new TGeoManager("geom", "FLUKA geometry");
197 fApplication->ConstructGeometry();
d59acfe7 198 if (!gGeoManager->IsClosed()) {
199 TGeoVolume *top = (TGeoVolume*)gGeoManager->GetListOfVolumes()->First();
200 gGeoManager->SetTopVolume(top);
201 gGeoManager->CloseGeometry("di");
202 } else {
203 TGeoNodeCache *cache = gGeoManager->GetCache();
204 if (!cache->HasIdArray()) {
a9ea1616 205 Warning("Init", "Node ID tracking must be enabled with TFluka: enabling...\n");
d59acfe7 206 cache->BuildIdArray();
207 }
208 }
829fb838 209 fNVolumes = fGeom->NofVolumes();
210 fGeom->CreateFlukaMatFile("flukaMat.inp");
211 if (fVerbosityLevel >=3) {
212 printf("== Number of volumes: %i\n ==", fNVolumes);
213 cout << "\t* InitPhysics() - Prepare input file to be called" << endl;
6d184c54 214 }
881cb248 215
216 fApplication->InitGeometry();
217
829fb838 218}
219
220
221//______________________________________________________________________________
222void TFluka::FinishGeometry() {
223//
224// Build-up table with region to medium correspondance
225//
226 if (fVerbosityLevel >=3) {
227 cout << "==> TFluka::FinishGeometry() called." << endl;
228 printf("----FinishGeometry - nothing to do with TGeo\n");
229 cout << "<== TFluka::FinishGeometry() called." << endl;
230 }
231}
232
233//______________________________________________________________________________
234void TFluka::BuildPhysics() {
235//
236// Prepare FLUKA input files and call FLUKA physics initialisation
237//
238
239 if (fVerbosityLevel >=3)
240 cout << "==> TFluka::BuildPhysics() called." << endl;
6d184c54 241
242
243 if (fVerbosityLevel >=3) {
244 TList *medlist = gGeoManager->GetListOfMedia();
245 TIter next(medlist);
246 TGeoMedium* med = 0x0;
247 TGeoMaterial* mat = 0x0;
248 Int_t ic = 0;
249
250 while((med = (TGeoMedium*)next()))
251 {
252 mat = med->GetMaterial();
253 printf("Medium %5d %12s %5d %5d\n", ic, (med->GetName()), med->GetId(), mat->GetIndex());
254 ic++;
255 }
256 }
257
258 //
259 // At this stage we have the information on materials and cuts available.
260 // Now create the pemf file
261
262 if (fGeneratePemf) fGeom->CreatePemfFile();
263
264 //
265 // Prepare input file with the current physics settings
266
829fb838 267 InitPhysics();
6d184c54 268
829fb838 269 cout << "\t* InitPhysics() - Prepare input file was called" << endl;
270
271 if (fVerbosityLevel >=2)
272 cout << "\t* Changing lfdrtr = (" << (GLOBAL.lfdrtr?'T':'F')
273 << ") in fluka..." << endl;
274 GLOBAL.lfdrtr = true;
275
276 if (fVerbosityLevel >=2)
277 cout << "\t* Opening file " << fInputFileName << endl;
278 const char* fname = fInputFileName;
6d184c54 279
829fb838 280 fluka_openinp(lunin, PASSCHARA(fname));
8e5bf079 281 fluka_openout(11, PASSCHARA("fluka.out"));
829fb838 282
283 if (fVerbosityLevel >=2)
284 cout << "\t* Calling flukam..." << endl;
285 flukam(1);
286
287 if (fVerbosityLevel >=2)
288 cout << "\t* Closing file " << fInputFileName << endl;
289 fluka_closeinp(lunin);
290
291 FinishGeometry();
292
293 if (fVerbosityLevel >=3)
294 cout << "<== TFluka::Init() called." << endl;
295
829fb838 296 if (fVerbosityLevel >=3)
297 cout << "<== TFluka::BuildPhysics() called." << endl;
298}
299
300//______________________________________________________________________________
301void TFluka::ProcessEvent() {
302//
303// Process one event
304//
b496f27c 305 if (fStopRun) {
a9ea1616 306 Warning("ProcessEvent", "User Run Abortion: No more events handled !\n");
b496f27c 307 fNEvent += 1;
308 return;
309 }
310
311 if (fVerbosityLevel >=3)
312 cout << "==> TFluka::ProcessEvent() called." << endl;
313 fApplication->GeneratePrimaries();
81f1d030 314 SOURCM.lsouit = true;
b496f27c 315 flukam(1);
316 if (fVerbosityLevel >=3)
317 cout << "<== TFluka::ProcessEvent() called." << endl;
318 //
319 // Increase event number
320 //
321 fNEvent += 1;
829fb838 322}
323
324//______________________________________________________________________________
325Bool_t TFluka::ProcessRun(Int_t nevent) {
326//
327// Run steering
328//
329
330 if (fVerbosityLevel >=3)
331 cout << "==> TFluka::ProcessRun(" << nevent << ") called."
332 << endl;
333
334 if (fVerbosityLevel >=2) {
335 cout << "\t* GLOBAL.fdrtr = " << (GLOBAL.lfdrtr?'T':'F') << endl;
336 cout << "\t* Calling flukam again..." << endl;
337 }
338
829fb838 339 Int_t todo = TMath::Abs(nevent);
340 for (Int_t ev = 0; ev < todo; ev++) {
341 fApplication->BeginEvent();
342 ProcessEvent();
343 fApplication->FinishEvent();
344 }
345
346 if (fVerbosityLevel >=3)
347 cout << "<== TFluka::ProcessRun(" << nevent << ") called."
348 << endl;
eea53470 349 // Write fluka specific scoring output
350 newplo();
351
829fb838 352 return kTRUE;
353}
354
355//_____________________________________________________________________________
356// methods for building/management of geometry
357
358// functions from GCONS
359//____________________________________________________________________________
360void TFluka::Gfmate(Int_t imat, char *name, Float_t &a, Float_t &z,
361 Float_t &dens, Float_t &radl, Float_t &absl,
362 Float_t* /*ubuf*/, Int_t& /*nbuf*/) {
363//
364 TGeoMaterial *mat;
365 TIter next (gGeoManager->GetListOfMaterials());
366 while ((mat = (TGeoMaterial*)next())) {
367 if (mat->GetUniqueID() == (UInt_t)imat) break;
368 }
369 if (!mat) {
370 Error("Gfmate", "no material with index %i found", imat);
371 return;
372 }
373 sprintf(name, "%s", mat->GetName());
374 a = mat->GetA();
375 z = mat->GetZ();
376 dens = mat->GetDensity();
377 radl = mat->GetRadLen();
378 absl = mat->GetIntLen();
379}
380
381//______________________________________________________________________________
382void TFluka::Gfmate(Int_t imat, char *name, Double_t &a, Double_t &z,
383 Double_t &dens, Double_t &radl, Double_t &absl,
384 Double_t* /*ubuf*/, Int_t& /*nbuf*/) {
385//
386 TGeoMaterial *mat;
387 TIter next (gGeoManager->GetListOfMaterials());
388 while ((mat = (TGeoMaterial*)next())) {
389 if (mat->GetUniqueID() == (UInt_t)imat) break;
390 }
391 if (!mat) {
392 Error("Gfmate", "no material with index %i found", imat);
393 return;
394 }
395 sprintf(name, "%s", mat->GetName());
396 a = mat->GetA();
397 z = mat->GetZ();
398 dens = mat->GetDensity();
399 radl = mat->GetRadLen();
400 absl = mat->GetIntLen();
401}
402
403// detector composition
404//______________________________________________________________________________
405void TFluka::Material(Int_t& kmat, const char* name, Double_t a,
406 Double_t z, Double_t dens, Double_t radl, Double_t absl,
407 Float_t* buf, Int_t nwbuf) {
408//
409 Double_t* dbuf = fGeom->CreateDoubleArray(buf, nwbuf);
410 Material(kmat, name, a, z, dens, radl, absl, dbuf, nwbuf);
411 delete [] dbuf;
412}
413
414//______________________________________________________________________________
415void TFluka::Material(Int_t& kmat, const char* name, Double_t a,
416 Double_t z, Double_t dens, Double_t radl, Double_t absl,
417 Double_t* /*buf*/, Int_t /*nwbuf*/) {
418//
fb2cbbec 419// Define a material
829fb838 420 TGeoMaterial *mat;
421 kmat = gGeoManager->GetListOfMaterials()->GetSize();
422 if ((z-Int_t(z)) > 1E-3) {
423 mat = fGeom->GetMakeWrongMaterial(z);
424 if (mat) {
425 mat->SetRadLen(radl,absl);
426 mat->SetUniqueID(kmat);
427 return;
428 }
429 }
430 gGeoManager->Material(name, a, z, dens, kmat, radl, absl);
431}
432
433//______________________________________________________________________________
434void TFluka::Mixture(Int_t& kmat, const char *name, Float_t *a,
435 Float_t *z, Double_t dens, Int_t nlmat, Float_t *wmat) {
fb2cbbec 436//
437// Define a material mixture
829fb838 438//
439 Double_t* da = fGeom->CreateDoubleArray(a, TMath::Abs(nlmat));
440 Double_t* dz = fGeom->CreateDoubleArray(z, TMath::Abs(nlmat));
441 Double_t* dwmat = fGeom->CreateDoubleArray(wmat, TMath::Abs(nlmat));
442
443 Mixture(kmat, name, da, dz, dens, nlmat, dwmat);
444 for (Int_t i=0; i<nlmat; i++) {
445 a[i] = da[i]; z[i] = dz[i]; wmat[i] = dwmat[i];
446 }
447
448 delete [] da;
449 delete [] dz;
450 delete [] dwmat;
451}
452
453//______________________________________________________________________________
454void TFluka::Mixture(Int_t& kmat, const char *name, Double_t *a,
455 Double_t *z, Double_t dens, Int_t nlmat, Double_t *wmat) {
456//
457 // Defines mixture OR COMPOUND IMAT as composed by
458 // THE BASIC NLMAT materials defined by arrays A,Z and WMAT
459 //
460 // If NLMAT > 0 then wmat contains the proportion by
461 // weights of each basic material in the mixture.
462 //
463 // If nlmat < 0 then WMAT contains the number of atoms
464 // of a given kind into the molecule of the COMPOUND
465 // In this case, WMAT in output is changed to relative
466 // weigths.
467 //
468 Int_t i,j;
469 if (nlmat < 0) {
470 nlmat = - nlmat;
471 Double_t amol = 0;
472 for (i=0;i<nlmat;i++) {
473 amol += a[i]*wmat[i];
474 }
475 for (i=0;i<nlmat;i++) {
476 wmat[i] *= a[i]/amol;
477 }
478 }
479 kmat = gGeoManager->GetListOfMaterials()->GetSize();
480 // Check if we have elements with fractional Z
481 TGeoMaterial *mat = 0;
482 TGeoMixture *mix = 0;
483 Bool_t mixnew = kFALSE;
484 for (i=0; i<nlmat; i++) {
485 if (z[i]-Int_t(z[i]) < 1E-3) continue;
486 // We have found an element with fractional Z -> loop mixtures to look for it
487 for (j=0; j<kmat; j++) {
488 mat = (TGeoMaterial*)gGeoManager->GetListOfMaterials()->At(j);
489 if (!mat) break;
490 if (!mat->IsMixture()) continue;
491 mix = (TGeoMixture*)mat;
492 if (TMath::Abs(z[i]-mix->GetZ()) >1E-3) continue;
829fb838 493 mixnew = kTRUE;
494 break;
495 }
496 if (!mixnew) Warning("Mixture","%s : cannot find component %i with fractional Z=%f\n", name, i, z[i]);
497 break;
498 }
499 if (mixnew) {
500 Int_t nlmatnew = nlmat+mix->GetNelements()-1;
501 Double_t *anew = new Double_t[nlmatnew];
502 Double_t *znew = new Double_t[nlmatnew];
503 Double_t *wmatnew = new Double_t[nlmatnew];
504 Int_t ind=0;
505 for (j=0; j<nlmat; j++) {
506 if (j==i) continue;
507 anew[ind] = a[j];
508 znew[ind] = z[j];
509 wmatnew[ind] = wmat[j];
510 ind++;
511 }
512 for (j=0; j<mix->GetNelements(); j++) {
513 anew[ind] = mix->GetAmixt()[j];
514 znew[ind] = mix->GetZmixt()[j];
515 wmatnew[ind] = wmat[i]*mix->GetWmixt()[j];
516 ind++;
517 }
518 Mixture(kmat, name, anew, znew, dens, nlmatnew, wmatnew);
519 delete [] anew;
520 delete [] znew;
521 delete [] wmatnew;
522 return;
523 }
524 // Now we need to compact identical elements within the mixture
525 // First check if this happens
526 mixnew = kFALSE;
527 for (i=0; i<nlmat-1; i++) {
528 for (j=i+1; j<nlmat; j++) {
529 if (z[i] == z[j]) {
530 mixnew = kTRUE;
531 break;
532 }
533 }
534 if (mixnew) break;
535 }
536 if (mixnew) {
537 Int_t nlmatnew = 0;
538 Double_t *anew = new Double_t[nlmat];
539 Double_t *znew = new Double_t[nlmat];
540 memset(znew, 0, nlmat*sizeof(Double_t));
541 Double_t *wmatnew = new Double_t[nlmat];
542 Bool_t skipi;
543 for (i=0; i<nlmat; i++) {
544 skipi = kFALSE;
545 for (j=0; j<nlmatnew; j++) {
546 if (z[i] == z[j]) {
547 wmatnew[j] += wmat[i];
548 skipi = kTRUE;
549 break;
550 }
551 }
552 if (skipi) continue;
553 anew[nlmatnew] = a[i];
554 znew[nlmatnew] = z[i];
555 wmatnew[nlmatnew] = wmat[i];
556 nlmatnew++;
557 }
558 Mixture(kmat, name, anew, znew, dens, nlmatnew, wmatnew);
559 delete [] anew;
560 delete [] znew;
561 delete [] wmatnew;
562 return;
563 }
564 gGeoManager->Mixture(name, a, z, dens, nlmat, wmat, kmat);
565}
566
567//______________________________________________________________________________
568void TFluka::Medium(Int_t& kmed, const char *name, Int_t nmat,
569 Int_t isvol, Int_t ifield, Double_t fieldm, Double_t tmaxfd,
570 Double_t stemax, Double_t deemax, Double_t epsil,
571 Double_t stmin, Float_t* ubuf, Int_t nbuf) {
b2129742 572 // Define a medium
573 //
829fb838 574 kmed = gGeoManager->GetListOfMedia()->GetSize()+1;
575 fMCGeo->Medium(kmed, name, nmat, isvol, ifield, fieldm, tmaxfd, stemax, deemax,
576 epsil, stmin, ubuf, nbuf);
577}
578
579//______________________________________________________________________________
580void TFluka::Medium(Int_t& kmed, const char *name, Int_t nmat,
581 Int_t isvol, Int_t ifield, Double_t fieldm, Double_t tmaxfd,
582 Double_t stemax, Double_t deemax, Double_t epsil,
583 Double_t stmin, Double_t* ubuf, Int_t nbuf) {
b2129742 584 // Define a medium
585 //
829fb838 586 kmed = gGeoManager->GetListOfMedia()->GetSize()+1;
587 fMCGeo->Medium(kmed, name, nmat, isvol, ifield, fieldm, tmaxfd, stemax, deemax,
588 epsil, stmin, ubuf, nbuf);
589}
590
591//______________________________________________________________________________
592void TFluka::Matrix(Int_t& krot, Double_t thetaX, Double_t phiX,
593 Double_t thetaY, Double_t phiY, Double_t thetaZ,
594 Double_t phiZ) {
595//
596 krot = gGeoManager->GetListOfMatrices()->GetEntriesFast();
597 fMCGeo->Matrix(krot, thetaX, phiX, thetaY, phiY, thetaZ, phiZ);
598}
599
600//______________________________________________________________________________
601void TFluka::Gstpar(Int_t itmed, const char* param, Double_t parval) {
602//
603//
7b203b6e 604//
829fb838 605 Bool_t process = kFALSE;
606 if (strncmp(param, "DCAY", 4) == 0 ||
607 strncmp(param, "PAIR", 4) == 0 ||
608 strncmp(param, "COMP", 4) == 0 ||
609 strncmp(param, "PHOT", 4) == 0 ||
610 strncmp(param, "PFIS", 4) == 0 ||
611 strncmp(param, "DRAY", 4) == 0 ||
612 strncmp(param, "ANNI", 4) == 0 ||
613 strncmp(param, "BREM", 4) == 0 ||
614 strncmp(param, "MUNU", 4) == 0 ||
615 strncmp(param, "CKOV", 4) == 0 ||
616 strncmp(param, "HADR", 4) == 0 ||
617 strncmp(param, "LOSS", 4) == 0 ||
618 strncmp(param, "MULS", 4) == 0 ||
619 strncmp(param, "RAYL", 4) == 0)
620 {
621 process = kTRUE;
622 }
81f1d030 623
829fb838 624 if (process) {
81f1d030 625 SetProcess(param, Int_t (parval), itmed);
829fb838 626 } else {
81f1d030 627 SetCut(param, parval, itmed);
829fb838 628 }
629}
630
631// functions from GGEOM
632//_____________________________________________________________________________
633void TFluka::Gsatt(const char *name, const char *att, Int_t val)
634{
6f5667d1 635 // Set visualisation attributes for one volume
829fb838 636 char vname[5];
637 fGeom->Vname(name,vname);
638 char vatt[5];
639 fGeom->Vname(att,vatt);
640 gGeoManager->SetVolumeAttribute(vname, vatt, val);
641}
642
643//______________________________________________________________________________
644Int_t TFluka::Gsvolu(const char *name, const char *shape, Int_t nmed,
645 Float_t *upar, Int_t np) {
646//
647 return fMCGeo->Gsvolu(name, shape, nmed, upar, np);
648}
649
650//______________________________________________________________________________
651Int_t TFluka::Gsvolu(const char *name, const char *shape, Int_t nmed,
652 Double_t *upar, Int_t np) {
653//
654 return fMCGeo->Gsvolu(name, shape, nmed, upar, np);
655}
656
657//______________________________________________________________________________
658void TFluka::Gsdvn(const char *name, const char *mother, Int_t ndiv,
659 Int_t iaxis) {
660//
661 fMCGeo->Gsdvn(name, mother, ndiv, iaxis);
662}
663
664//______________________________________________________________________________
665void TFluka::Gsdvn2(const char *name, const char *mother, Int_t ndiv,
666 Int_t iaxis, Double_t c0i, Int_t numed) {
667//
668 fMCGeo->Gsdvn2(name, mother, ndiv, iaxis, c0i, numed);
669}
670
671//______________________________________________________________________________
672void TFluka::Gsdvt(const char *name, const char *mother, Double_t step,
673 Int_t iaxis, Int_t numed, Int_t ndvmx) {
674//
675 fMCGeo->Gsdvt(name, mother, step, iaxis, numed, ndvmx);
676}
677
678//______________________________________________________________________________
679void TFluka::Gsdvt2(const char *name, const char *mother, Double_t step,
680 Int_t iaxis, Double_t c0, Int_t numed, Int_t ndvmx) {
681//
682 fMCGeo->Gsdvt2(name, mother, step, iaxis, c0, numed, ndvmx);
683}
684
685//______________________________________________________________________________
686void TFluka::Gsord(const char * /*name*/, Int_t /*iax*/) {
687//
688// Nothing to do with TGeo
689}
690
691//______________________________________________________________________________
692void TFluka::Gspos(const char *name, Int_t nr, const char *mother,
693 Double_t x, Double_t y, Double_t z, Int_t irot,
694 const char *konly) {
695//
696 fMCGeo->Gspos(name, nr, mother, x, y, z, irot, konly);
697}
698
699//______________________________________________________________________________
700void TFluka::Gsposp(const char *name, Int_t nr, const char *mother,
701 Double_t x, Double_t y, Double_t z, Int_t irot,
702 const char *konly, Float_t *upar, Int_t np) {
703 //
704 fMCGeo->Gsposp(name, nr, mother, x, y, z, irot, konly, upar, np);
705}
706
707//______________________________________________________________________________
708void TFluka::Gsposp(const char *name, Int_t nr, const char *mother,
709 Double_t x, Double_t y, Double_t z, Int_t irot,
710 const char *konly, Double_t *upar, Int_t np) {
711 //
712 fMCGeo->Gsposp(name, nr, mother, x, y, z, irot, konly, upar, np);
713}
714
715//______________________________________________________________________________
716void TFluka::Gsbool(const char* /*onlyVolName*/, const char* /*manyVolName*/) {
717//
718// Nothing to do with TGeo
719}
720
a9ea1616 721//______________________________________________________________________
722Bool_t TFluka::GetTransformation(const TString &volumePath,TGeoHMatrix &mat)
723{
724 // Returns the Transformation matrix between the volume specified
725 // by the path volumePath and the Top or mater volume. The format
726 // of the path volumePath is as follows (assuming ALIC is the Top volume)
727 // "/ALIC_1/DDIP_1/S05I_2/S05H_1/S05G_3". Here ALIC is the top most
728 // or master volume which has only 1 instance of. Of all of the daughter
729 // volumes of ALICE, DDIP volume copy #1 is indicated. Similarly for
730 // the daughter volume of DDIP is S05I copy #2 and so on.
731 // Inputs:
732 // TString& volumePath The volume path to the specific volume
733 // for which you want the matrix. Volume name
734 // hierarchy is separated by "/" while the
735 // copy number is appended using a "_".
736 // Outputs:
737 // TGeoHMatrix &mat A matrix with its values set to those
738 // appropriate to the Local to Master transformation
739 // Return:
740 // A logical value if kFALSE then an error occurred and no change to
741 // mat was made.
742
743 // We have to preserve the modeler state
744 return fMCGeo->GetTransformation(volumePath, mat);
745}
746
747//______________________________________________________________________
748Bool_t TFluka::GetShape(const TString &volumePath,TString &shapeType,
749 TArrayD &par)
750{
751 // Returns the shape and its parameters for the volume specified
752 // by volumeName.
753 // Inputs:
754 // TString& volumeName The volume name
755 // Outputs:
756 // TString &shapeType Shape type
757 // TArrayD &par A TArrayD of parameters with all of the
758 // parameters of the specified shape.
759 // Return:
760 // A logical indicating whether there was an error in getting this
761 // information
762 return fMCGeo->GetShape(volumePath, shapeType, par);
763}
764
765//______________________________________________________________________
766Bool_t TFluka::GetMaterial(const TString &volumeName,
767 TString &name,Int_t &imat,
768 Double_t &a,Double_t &z,Double_t &dens,
769 Double_t &radl,Double_t &inter,TArrayD &par)
770{
771 // Returns the Material and its parameters for the volume specified
772 // by volumeName.
773 // Note, Geant3 stores and uses mixtures as an element with an effective
774 // Z and A. Consequently, if the parameter Z is not integer, then
775 // this material represents some sort of mixture.
776 // Inputs:
777 // TString& volumeName The volume name
778 // Outputs:
779 // TSrting &name Material name
780 // Int_t &imat Material index number
781 // Double_t &a Average Atomic mass of material
782 // Double_t &z Average Atomic number of material
783 // Double_t &dens Density of material [g/cm^3]
784 // Double_t &radl Average radiation length of material [cm]
785 // Double_t &inter Average interaction length of material [cm]
786 // TArrayD &par A TArrayD of user defined parameters.
787 // Return:
788 // kTRUE if no errors
789 return fMCGeo->GetMaterial(volumeName,name,imat,a,z,dens,radl,inter,par);
790}
791
792//______________________________________________________________________
793Bool_t TFluka::GetMedium(const TString &volumeName,TString &name,
794 Int_t &imed,Int_t &nmat,Int_t &isvol,Int_t &ifield,
795 Double_t &fieldm,Double_t &tmaxfd,Double_t &stemax,
796 Double_t &deemax,Double_t &epsil, Double_t &stmin,
797 TArrayD &par)
798{
799 // Returns the Medium and its parameters for the volume specified
800 // by volumeName.
801 // Inputs:
802 // TString& volumeName The volume name.
803 // Outputs:
804 // TString &name Medium name
805 // Int_t &nmat Material number defined for this medium
806 // Int_t &imed The medium index number
807 // Int_t &isvol volume number defined for this medium
808 // Int_t &iflield Magnetic field flag
809 // Double_t &fieldm Magnetic field strength
810 // Double_t &tmaxfd Maximum angle of deflection per step
811 // Double_t &stemax Maximum step size
812 // Double_t &deemax Maximum fraction of energy allowed to be lost
813 // to continuous process.
814 // Double_t &epsil Boundary crossing precision
815 // Double_t &stmin Minimum step size allowed
816 // TArrayD &par A TArrayD of user parameters with all of the
817 // parameters of the specified medium.
818 // Return:
819 // kTRUE if there where no errors
820 return fMCGeo->GetMedium(volumeName,name,imed,nmat,isvol,ifield,fieldm,tmaxfd,stemax,deemax,epsil,stmin,par);
821}
822
829fb838 823//______________________________________________________________________________
824void TFluka::SetCerenkov(Int_t itmed, Int_t npckov, Float_t* ppckov,
825 Float_t* absco, Float_t* effic, Float_t* rindex) {
826//
827// Set Cerenkov properties for medium itmed
828//
829// npckov: number of sampling points
830// ppckov: energy values
831// absco: absorption length
832// effic: quantum efficiency
833// rindex: refraction index
834//
835//
836//
837// Create object holding Cerenkov properties
838//
839 TFlukaCerenkov* cerenkovProperties = new TFlukaCerenkov(npckov, ppckov, absco, effic, rindex);
840//
841// Pass object to medium
842 TGeoMedium* medium = gGeoManager->GetMedium(itmed);
843 medium->SetCerenkovProperties(cerenkovProperties);
844}
845
b2be0e73 846void TFluka::SetCerenkov(Int_t itmed, Int_t npckov, Float_t* ppckov,
847 Float_t* absco, Float_t* effic, Float_t* rindex, Float_t* rfl) {
848//
849// Set Cerenkov properties for medium itmed
850//
851// npckov: number of sampling points
852// ppckov: energy values
853// absco: absorption length
854// effic: quantum efficiency
855// rindex: refraction index
856// rfl: reflectivity for boundary to medium itmed
857//
858//
859// Create object holding Cerenkov properties
860//
861 TFlukaCerenkov* cerenkovProperties = new TFlukaCerenkov(npckov, ppckov, absco, effic, rindex, rfl);
862//
863// Pass object to medium
864 TGeoMedium* medium = gGeoManager->GetMedium(itmed);
865 medium->SetCerenkovProperties(cerenkovProperties);
866}
867
868
829fb838 869//______________________________________________________________________________
870void TFluka::SetCerenkov(Int_t /*itmed*/, Int_t /*npckov*/, Double_t * /*ppckov*/,
871 Double_t * /*absco*/, Double_t * /*effic*/, Double_t * /*rindex*/) {
872//
b2be0e73 873// Double_t version not implemented
829fb838 874}
b2be0e73 875
876void TFluka::SetCerenkov(Int_t /*itmed*/, Int_t /*npckov*/, Double_t* /*ppckov*/,
877 Double_t* /*absco*/, Double_t* /*effic*/, Double_t* /*rindex*/, Double_t* /*rfl*/) {
878//
879// // Double_t version not implemented
880}
881
829fb838 882// Euclid
883//______________________________________________________________________________
884void TFluka::WriteEuclid(const char* /*fileName*/, const char* /*topVol*/,
885 Int_t /*number*/, Int_t /*nlevel*/) {
886//
887// Not with TGeo
a9ea1616 888 Warning("WriteEuclid", "Not implemented !");
829fb838 889}
890
891
892
893//_____________________________________________________________________________
894// methods needed by the stepping
895//____________________________________________________________________________
896
897Int_t TFluka::GetMedium() const {
898//
899// Get the medium number for the current fluka region
900//
901 return fGeom->GetMedium(); // this I need to check due to remapping !!!
902}
903
a9ea1616 904//____________________________________________________________________________
905Int_t TFluka::GetDummyRegion() const
906{
907// Returns index of the dummy region.
908 return fGeom->GetDummyRegion();
909}
829fb838 910
a9ea1616 911//____________________________________________________________________________
912Int_t TFluka::GetDummyLattice() const
913{
914// Returns index of the dummy lattice.
915 return fGeom->GetDummyLattice();
916}
829fb838 917
918//____________________________________________________________________________
919// particle table usage
920// ID <--> PDG transformations
921//_____________________________________________________________________________
922Int_t TFluka::IdFromPDG(Int_t pdg) const
923{
924 //
925 // Return Fluka code from PDG and pseudo ENDF code
926
927 // Catch the feedback photons
a9ea1616 928 if (pdg == 50000051) return (kFLUKAoptical);
829fb838 929 // MCIHAD() goes from pdg to fluka internal.
930 Int_t intfluka = mcihad(pdg);
931 // KPTOIP array goes from internal to official
932 return GetFlukaKPTOIP(intfluka);
933}
934
935//______________________________________________________________________________
936Int_t TFluka::PDGFromId(Int_t id) const
937{
938 //
939 // Return PDG code and pseudo ENDF code from Fluka code
f926898e 940 // Alpha He3 Triton Deuteron gen. ion opt. photon
941 Int_t idSpecial[6] = {10020040, 10020030, 10010030, 10010020, 10000000, 50000050};
829fb838 942 // IPTOKP array goes from official to internal
943
a9ea1616 944 if (id == kFLUKAoptical) {
829fb838 945// Cerenkov photon
bd3d5c8a 946 if (fVerbosityLevel >= 3)
829fb838 947 printf("\n PDGFromId: Cerenkov Photon \n");
948 return 50000050;
949 }
950// Error id
ece92b30 951 if (id == 0 || id < kFLUKAcodemin || id > kFLUKAcodemax) {
f926898e 952 if (fVerbosityLevel >= 3)
829fb838 953 printf("PDGFromId: Error id = 0\n");
954 return -1;
955 }
956// Good id
f926898e 957 if (id > 0) {
958 Int_t intfluka = GetFlukaIPTOKP(id);
959 if (intfluka == 0) {
960 if (fVerbosityLevel >= 3)
961 printf("PDGFromId: Error intfluka = 0: %d\n", id);
962 return -1;
963 } else if (intfluka < 0) {
964 if (fVerbosityLevel >= 3)
965 printf("PDGFromId: Error intfluka < 0: %d\n", id);
966 return -1;
967 }
968 if (fVerbosityLevel >= 3)
969 printf("mpdgha called with %d %d \n", id, intfluka);
f926898e 970 return mpdgha(intfluka);
971 } else {
972 // ions and optical photons
ece92b30 973 return idSpecial[id - kFLUKAcodemin];
829fb838 974 }
829fb838 975}
976
bd3d5c8a 977void TFluka::StopTrack()
978{
979 // Set stopping conditions
980 // Works for photons and charged particles
981 fStopped = kTRUE;
982}
983
829fb838 984//_____________________________________________________________________________
985// methods for physics management
986//____________________________________________________________________________
987//
988// set methods
989//
990
1df5fa54 991void TFluka::SetProcess(const char* flagName, Int_t flagValue, Int_t imed)
829fb838 992{
993// Set process user flag for material imat
1df5fa54 994//
995//
996// Update if already in the list
829fb838 997//
fb2cbbec 998 TIter next(fUserConfig);
1df5fa54 999 TFlukaConfigOption* proc;
1000 while((proc = (TFlukaConfigOption*)next()))
1001 {
fb2cbbec 1002 if (proc->Medium() == imed) {
1003 proc->SetProcess(flagName, flagValue);
1004 return;
1005 }
1df5fa54 1006 }
fb2cbbec 1007 proc = new TFlukaConfigOption(imed);
1008 proc->SetProcess(flagName, flagValue);
1009 fUserConfig->Add(proc);
1010}
1011
1012//______________________________________________________________________________
1013Bool_t TFluka::SetProcess(const char* flagName, Int_t flagValue)
1014{
1015// Set process user flag
1df5fa54 1016//
1df5fa54 1017//
fb2cbbec 1018 SetProcess(flagName, flagValue, -1);
1df5fa54 1019 return kTRUE;
829fb838 1020}
1021
1022//______________________________________________________________________________
1023void TFluka::SetCut(const char* cutName, Double_t cutValue, Int_t imed)
1024{
1025// Set user cut value for material imed
1026//
fb2cbbec 1027 TIter next(fUserConfig);
1028 TFlukaConfigOption* proc;
1029 while((proc = (TFlukaConfigOption*)next()))
1030 {
1031 if (proc->Medium() == imed) {
1032 proc->SetCut(cutName, cutValue);
1033 return;
1034 }
1035 }
1036
1037 proc = new TFlukaConfigOption(imed);
1038 proc->SetCut(cutName, cutValue);
1039 fUserConfig->Add(proc);
829fb838 1040}
1041
1042//______________________________________________________________________________
1043Bool_t TFluka::SetCut(const char* cutName, Double_t cutValue)
1044{
1045// Set user cut value
1046//
1df5fa54 1047//
fb2cbbec 1048 SetCut(cutName, cutValue, -1);
1049 return kTRUE;
829fb838 1050}
1051
f450e9d0 1052
1053void TFluka::SetUserScoring(const char* option, Int_t npr, char* outfile, Float_t* what)
b496f27c 1054{
1055//
f450e9d0 1056// Adds a user scoring option to the list
b496f27c 1057//
f450e9d0 1058 TFlukaScoringOption* opt = new TFlukaScoringOption(option, "User Scoring", npr,outfile,what);
1059 fUserScore->Add(opt);
1060}
1061//______________________________________________________________________________
1062void TFluka::SetUserScoring(const char* option, Int_t npr, char* outfile, Float_t* what, const char* det1, const char* det2, const char* det3)
1063{
1064//
1065// Adds a user scoring option to the list
1066//
1067 TFlukaScoringOption* opt = new TFlukaScoringOption(option, "User Scoring", npr, outfile, what, det1, det2, det3);
b496f27c 1068 fUserScore->Add(opt);
1069}
b496f27c 1070
829fb838 1071//______________________________________________________________________________
1072Double_t TFluka::Xsec(char*, Double_t, Int_t, Int_t)
1073{
a9ea1616 1074 Warning("Xsec", "Not yet implemented.!\n"); return -1.;
829fb838 1075}
1076
1077
1078//______________________________________________________________________________
1079void TFluka::InitPhysics()
1080{
1081//
1082// Physics initialisation with preparation of FLUKA input cards
1083//
fb2cbbec 1084// Construct file names
1085 FILE *pFlukaVmcCoreInp, *pFlukaVmcFlukaMat, *pFlukaVmcInp;
1086 TString sFlukaVmcCoreInp = getenv("ALICE_ROOT");
1087 sFlukaVmcCoreInp +="/TFluka/input/";
1088 TString sFlukaVmcTmp = "flukaMat.inp";
1089 TString sFlukaVmcInp = GetInputFileName();
1090 sFlukaVmcCoreInp += GetCoreInputFileName();
1091
1092// Open files
1093 if ((pFlukaVmcCoreInp = fopen(sFlukaVmcCoreInp.Data(),"r")) == NULL) {
a9ea1616 1094 Warning("InitPhysics", "\nCannot open file %s\n",sFlukaVmcCoreInp.Data());
fb2cbbec 1095 exit(1);
1096 }
1097 if ((pFlukaVmcFlukaMat = fopen(sFlukaVmcTmp.Data(),"r")) == NULL) {
a9ea1616 1098 Warning("InitPhysics", "\nCannot open file %s\n",sFlukaVmcTmp.Data());
fb2cbbec 1099 exit(1);
1100 }
1101 if ((pFlukaVmcInp = fopen(sFlukaVmcInp.Data(),"w")) == NULL) {
a9ea1616 1102 Warning("InitPhysics", "\nCannot open file %s\n",sFlukaVmcInp.Data());
fb2cbbec 1103 exit(1);
1104 }
829fb838 1105
fb2cbbec 1106// Copy core input file
1107 Char_t sLine[255];
1108 Float_t fEventsPerRun;
829fb838 1109
fb2cbbec 1110 while ((fgets(sLine,255,pFlukaVmcCoreInp)) != NULL) {
1111 if (strncmp(sLine,"GEOEND",6) != 0)
1112 fprintf(pFlukaVmcInp,"%s",sLine); // copy until GEOEND card
1113 else {
1114 fprintf(pFlukaVmcInp,"GEOEND\n"); // add GEOEND card
1115 goto flukamat;
829fb838 1116 }
fb2cbbec 1117 } // end of while until GEOEND card
1118
829fb838 1119
fb2cbbec 1120 flukamat:
1121 while ((fgets(sLine,255,pFlukaVmcFlukaMat)) != NULL) { // copy flukaMat.inp file
1122 fprintf(pFlukaVmcInp,"%s\n",sLine);
1123 }
1124
1125 while ((fgets(sLine,255,pFlukaVmcCoreInp)) != NULL) {
1126 if (strncmp(sLine,"START",5) != 0)
1127 fprintf(pFlukaVmcInp,"%s\n",sLine);
1128 else {
1129 sscanf(sLine+10,"%10f",&fEventsPerRun);
1130 goto fin;
1131 }
1132 } //end of while until START card
1133
1134 fin:
829fb838 1135
f450e9d0 1136
1137// Pass information to configuration objects
829fb838 1138
fb2cbbec 1139 Float_t fLastMaterial = fGeom->GetLastMaterialIndex();
1140 TFlukaConfigOption::SetStaticInfo(pFlukaVmcInp, 3, fLastMaterial, fGeom);
1141
1142 TIter next(fUserConfig);
1143 TFlukaConfigOption* proc;
f450e9d0 1144 while((proc = dynamic_cast<TFlukaConfigOption*> (next()))) proc->WriteFlukaInputCards();
1145//
1146// Process Fluka specific scoring options
1147//
1148 TFlukaScoringOption::SetStaticInfo(pFlukaVmcInp, fGeom);
1149 Float_t loginp = 49.0;
1150 Int_t inp = 0;
1151 Int_t nscore = fUserScore->GetEntries();
1152
a9ea1616 1153 TFlukaScoringOption *mopo = 0;
1154 TFlukaScoringOption *mopi = 0;
fb2cbbec 1155
f450e9d0 1156 for (Int_t isc = 0; isc < nscore; isc++)
1157 {
1158 mopo = dynamic_cast<TFlukaScoringOption*> (fUserScore->At(isc));
1159 char* fileName = mopo->GetFileName();
1160 Int_t size = strlen(fileName);
1161 Float_t lun = -1.;
1162//
1163// Check if new output file has to be opened
1164 for (Int_t isci = 0; isci < isc; isci++) {
a9d74780 1165
1166
1167 mopi = dynamic_cast<TFlukaScoringOption*> (fUserScore->At(isci));
f450e9d0 1168 if(strncmp(mopi->GetFileName(), fileName, size)==0) {
1169 //
1170 // No, the file already exists
1171 lun = mopi->GetLun();
1172 mopo->SetLun(lun);
1173 break;
1174 }
1175 } // inner loop
1176
1177 if (lun == -1.) {
1178 // Open new output file
1179 inp++;
1180 mopo->SetLun(loginp + inp);
1181 mopo->WriteOpenFlukaFile();
1182 }
1183 mopo->WriteFlukaInputCards();
1184 }
1185
829fb838 1186// Add START and STOP card
f450e9d0 1187 fprintf(pFlukaVmcInp,"START %10.1f\n",fEventsPerRun);
1188 fprintf(pFlukaVmcInp,"STOP \n");
829fb838 1189
1190
1191// Close files
3b8c325d 1192 fclose(pFlukaVmcCoreInp);
1193 fclose(pFlukaVmcFlukaMat);
1194 fclose(pFlukaVmcInp);
fb2cbbec 1195
1196
1197//
1198// Initialisation needed for Cerenkov photon production and transport
1199 TObjArray *matList = GetFlukaMaterials();
1200 Int_t nmaterial = matList->GetEntriesFast();
1201 fMaterials = new Int_t[nmaterial+3];
1202
1203 for (Int_t im = 0; im < nmaterial; im++)
1204 {
1205 TGeoMaterial* material = dynamic_cast<TGeoMaterial*> (matList->At(im));
1206 Int_t idmat = material->GetIndex();
1207 fMaterials[idmat] = im;
1208 }
829fb838 1209} // end of InitPhysics
1210
1211
1212//______________________________________________________________________________
07f5b33e 1213void TFluka::SetMaxStep(Double_t step)
829fb838 1214{
07f5b33e 1215// Set the maximum step size
1216 if (step > 1.e4) return;
1217
1218 Int_t mreg, latt;
1219 fGeom->GetCurrentRegion(mreg, latt);
9c0c08ce 1220 STEPSZ.stepmx[mreg - 1] = step;
829fb838 1221}
1222
2f09b80e 1223
1224Double_t TFluka::MaxStep() const
1225{
1226// Return the maximum for current medium
1227 Int_t mreg, latt;
1228 fGeom->GetCurrentRegion(mreg, latt);
1229 return (STEPSZ.stepmx[mreg - 1]);
1230}
1231
829fb838 1232//______________________________________________________________________________
1233void TFluka::SetMaxNStep(Int_t)
1234{
1235// SetMaxNStep is dummy procedure in TFluka !
1236 if (fVerbosityLevel >=3)
1237 cout << "SetMaxNStep is dummy procedure in TFluka !" << endl;
1238}
1239
1240//______________________________________________________________________________
1241void TFluka::SetUserDecay(Int_t)
1242{
1243// SetUserDecay is dummy procedure in TFluka !
1244 if (fVerbosityLevel >=3)
1245 cout << "SetUserDecay is dummy procedure in TFluka !" << endl;
1246}
1247
1248//
1249// dynamic properties
1250//
1251//______________________________________________________________________________
1252void TFluka::TrackPosition(TLorentzVector& position) const
1253{
1254// Return the current position in the master reference frame of the
1255// track being transported
1256// TRACKR.atrack = age of the particle
1257// TRACKR.xtrack = x-position of the last point
1258// TRACKR.ytrack = y-position of the last point
1259// TRACKR.ztrack = z-position of the last point
a9ea1616 1260 FlukaCallerCode_t caller = GetCaller();
1261 if (caller == kENDRAW || caller == kUSDRAW ||
1262 caller == kBXExiting || caller == kBXEntering ||
1263 caller == kUSTCKV) {
829fb838 1264 position.SetX(GetXsco());
1265 position.SetY(GetYsco());
1266 position.SetZ(GetZsco());
1267 position.SetT(TRACKR.atrack);
1268 }
a9ea1616 1269 else if (caller == kMGDRAW) {
829fb838 1270 position.SetX(TRACKR.xtrack[TRACKR.ntrack]);
1271 position.SetY(TRACKR.ytrack[TRACKR.ntrack]);
1272 position.SetZ(TRACKR.ztrack[TRACKR.ntrack]);
1273 position.SetT(TRACKR.atrack);
1274 }
a9ea1616 1275 else if (caller == kSODRAW) {
829fb838 1276 position.SetX(TRACKR.xtrack[TRACKR.ntrack]);
1277 position.SetY(TRACKR.ytrack[TRACKR.ntrack]);
1278 position.SetZ(TRACKR.ztrack[TRACKR.ntrack]);
1279 position.SetT(0);
a9ea1616 1280 } else if (caller == kMGResumedTrack) {
5d80a015 1281 position.SetX(TRACKR.spausr[0]);
1282 position.SetY(TRACKR.spausr[1]);
1283 position.SetZ(TRACKR.spausr[2]);
1284 position.SetT(TRACKR.spausr[3]);
829fb838 1285 }
1286 else
1287 Warning("TrackPosition","position not available");
1288}
1289
1290//______________________________________________________________________________
1291void TFluka::TrackPosition(Double_t& x, Double_t& y, Double_t& z) const
1292{
1293// Return the current position in the master reference frame of the
1294// track being transported
1295// TRACKR.atrack = age of the particle
1296// TRACKR.xtrack = x-position of the last point
1297// TRACKR.ytrack = y-position of the last point
1298// TRACKR.ztrack = z-position of the last point
a9ea1616 1299 FlukaCallerCode_t caller = GetCaller();
1300 if (caller == kENDRAW || caller == kUSDRAW ||
1301 caller == kBXExiting || caller == kBXEntering ||
1302 caller == kUSTCKV) {
829fb838 1303 x = GetXsco();
1304 y = GetYsco();
1305 z = GetZsco();
1306 }
a9ea1616 1307 else if (caller == kMGDRAW || caller == kSODRAW) {
829fb838 1308 x = TRACKR.xtrack[TRACKR.ntrack];
1309 y = TRACKR.ytrack[TRACKR.ntrack];
1310 z = TRACKR.ztrack[TRACKR.ntrack];
1311 }
a9ea1616 1312 else if (caller == kMGResumedTrack) {
5d80a015 1313 x = TRACKR.spausr[0];
1314 y = TRACKR.spausr[1];
1315 z = TRACKR.spausr[2];
1316 }
829fb838 1317 else
1318 Warning("TrackPosition","position not available");
1319}
1320
1321//______________________________________________________________________________
1322void TFluka::TrackMomentum(TLorentzVector& momentum) const
1323{
1324// Return the direction and the momentum (GeV/c) of the track
1325// currently being transported
1326// TRACKR.ptrack = momentum of the particle (not always defined, if
1327// < 0 must be obtained from etrack)
1328// TRACKR.cx,y,ztrck = direction cosines of the current particle
1329// TRACKR.etrack = total energy of the particle
1330// TRACKR.jtrack = identity number of the particle
1331// PAPROP.am[TRACKR.jtrack] = particle mass in gev
a9ea1616 1332 FlukaCallerCode_t caller = GetCaller();
1333 FlukaProcessCode_t icode = GetIcode();
1334
1335 if (caller != kEEDRAW && caller != kMGResumedTrack &&
1336 (caller != kENDRAW || (icode != kEMFSCOstopping1 && icode != kEMFSCOstopping2))) {
829fb838 1337 if (TRACKR.ptrack >= 0) {
1338 momentum.SetPx(TRACKR.ptrack*TRACKR.cxtrck);
1339 momentum.SetPy(TRACKR.ptrack*TRACKR.cytrck);
1340 momentum.SetPz(TRACKR.ptrack*TRACKR.cztrck);
1341 momentum.SetE(TRACKR.etrack);
1342 return;
1343 }
1344 else {
ece92b30 1345 Double_t p = sqrt(TRACKR.etrack * TRACKR.etrack - ParticleMassFPC(TRACKR.jtrack) * ParticleMassFPC(TRACKR.jtrack));
829fb838 1346 momentum.SetPx(p*TRACKR.cxtrck);
1347 momentum.SetPy(p*TRACKR.cytrck);
1348 momentum.SetPz(p*TRACKR.cztrck);
1349 momentum.SetE(TRACKR.etrack);
1350 return;
1351 }
a9ea1616 1352 } else if (caller == kMGResumedTrack) {
5d80a015 1353 momentum.SetPx(TRACKR.spausr[4]);
1354 momentum.SetPy(TRACKR.spausr[5]);
1355 momentum.SetPz(TRACKR.spausr[6]);
1356 momentum.SetE (TRACKR.spausr[7]);
1357 return;
a9ea1616 1358 } else if (caller == kENDRAW && (icode == kEMFSCOstopping1 || icode == kEMFSCOstopping2)) {
1359 momentum.SetPx(0.);
1360 momentum.SetPy(0.);
1361 momentum.SetPz(0.);
1362 momentum.SetE(TrackMass());
829fb838 1363 }
1364 else
1365 Warning("TrackMomentum","momentum not available");
1366}
1367
1368//______________________________________________________________________________
1369void TFluka::TrackMomentum(Double_t& px, Double_t& py, Double_t& pz, Double_t& e) const
1370{
1371// Return the direction and the momentum (GeV/c) of the track
1372// currently being transported
1373// TRACKR.ptrack = momentum of the particle (not always defined, if
1374// < 0 must be obtained from etrack)
1375// TRACKR.cx,y,ztrck = direction cosines of the current particle
1376// TRACKR.etrack = total energy of the particle
1377// TRACKR.jtrack = identity number of the particle
1378// PAPROP.am[TRACKR.jtrack] = particle mass in gev
a9ea1616 1379 FlukaCallerCode_t caller = GetCaller();
1380 FlukaProcessCode_t icode = GetIcode();
1381 if (caller != kEEDRAW && caller != kMGResumedTrack &&
1382 (caller != kENDRAW || (icode != kEMFSCOstopping1 && icode != kEMFSCOstopping2))) {
829fb838 1383 if (TRACKR.ptrack >= 0) {
1384 px = TRACKR.ptrack*TRACKR.cxtrck;
1385 py = TRACKR.ptrack*TRACKR.cytrck;
1386 pz = TRACKR.ptrack*TRACKR.cztrck;
a9ea1616 1387 e = TRACKR.etrack;
829fb838 1388 return;
1389 }
1390 else {
ece92b30 1391 Double_t p = sqrt(TRACKR.etrack * TRACKR.etrack - ParticleMassFPC(TRACKR.jtrack) * ParticleMassFPC(TRACKR.jtrack));
829fb838 1392 px = p*TRACKR.cxtrck;
1393 py = p*TRACKR.cytrck;
1394 pz = p*TRACKR.cztrck;
a9ea1616 1395 e = TRACKR.etrack;
829fb838 1396 return;
1397 }
a9ea1616 1398 } else if (caller == kMGResumedTrack) {
5d80a015 1399 px = TRACKR.spausr[4];
1400 py = TRACKR.spausr[5];
1401 pz = TRACKR.spausr[6];
1402 e = TRACKR.spausr[7];
0773d0ac 1403 return;
a9ea1616 1404 } else if (caller == kENDRAW && (icode == kEMFSCOstopping1 || icode == kEMFSCOstopping2)) {
1405 px = 0.;
1406 py = 0.;
1407 pz = 0.;
1408 e = TrackMass();
829fb838 1409 }
1410 else
1411 Warning("TrackMomentum","momentum not available");
1412}
1413
1414//______________________________________________________________________________
1415Double_t TFluka::TrackStep() const
1416{
1417// Return the length in centimeters of the current step
1418// TRACKR.ctrack = total curved path
a9ea1616 1419 FlukaCallerCode_t caller = GetCaller();
1420 if (caller == kBXEntering || caller == kBXExiting ||
1421 caller == kENDRAW || caller == kUSDRAW ||
1422 caller == kUSTCKV || caller == kMGResumedTrack)
829fb838 1423 return 0.0;
a9ea1616 1424 else if (caller == kMGDRAW)
829fb838 1425 return TRACKR.ctrack;
669cede4 1426 else {
1427 Warning("TrackStep", "track step not available");
1428 return 0.0;
1429 }
829fb838 1430}
1431
1432//______________________________________________________________________________
1433Double_t TFluka::TrackLength() const
1434{
1435// TRACKR.cmtrck = cumulative curved path since particle birth
a9ea1616 1436 FlukaCallerCode_t caller = GetCaller();
1437 if (caller == kBXEntering || caller == kBXExiting ||
1438 caller == kENDRAW || caller == kUSDRAW || caller == kMGDRAW ||
1439 caller == kUSTCKV)
829fb838 1440 return TRACKR.cmtrck;
a9ea1616 1441 else if (caller == kMGResumedTrack)
5d80a015 1442 return TRACKR.spausr[8];
669cede4 1443 else {
1444 Warning("TrackLength", "track length not available");
1445 return 0.0;
1446 }
829fb838 1447}
1448
1449//______________________________________________________________________________
1450Double_t TFluka::TrackTime() const
1451{
1452// Return the current time of flight of the track being transported
1453// TRACKR.atrack = age of the particle
a9ea1616 1454 FlukaCallerCode_t caller = GetCaller();
1455 if (caller == kBXEntering || caller == kBXExiting ||
1456 caller == kENDRAW || caller == kUSDRAW || caller == kMGDRAW ||
1457 caller == kUSTCKV)
829fb838 1458 return TRACKR.atrack;
a9ea1616 1459 else if (caller == kMGResumedTrack)
5d80a015 1460 return TRACKR.spausr[3];
669cede4 1461 else {
1462 Warning("TrackTime", "track time not available");
1463 return 0.0;
1464 }
829fb838 1465}
1466
1467//______________________________________________________________________________
1468Double_t TFluka::Edep() const
1469{
1470// Energy deposition
1471// if TRACKR.ntrack = 0, TRACKR.mtrack = 0:
1472// -->local energy deposition (the value and the point are not recorded in TRACKR)
1473// but in the variable "rull" of the procedure "endraw.cxx"
1474// if TRACKR.ntrack > 0, TRACKR.mtrack = 0:
1475// -->no energy loss along the track
1476// if TRACKR.ntrack > 0, TRACKR.mtrack > 0:
1477// -->energy loss distributed along the track
07f5b33e 1478// TRACKR.dtrack = energy deposition of the jth deposition event
829fb838 1479
1480 // If coming from bxdraw we have 2 steps of 0 length and 0 edep
669cede4 1481 // If coming from usdraw we just signal particle production - no edep
1482 // If just first time after resuming, no edep for the primary
a9ea1616 1483 FlukaCallerCode_t caller = GetCaller();
1484 if (caller == kBXExiting || caller == kBXEntering ||
1485 caller == kUSDRAW || caller == kMGResumedTrack) return 0.0;
829fb838 1486 Double_t sum = 0;
1487 for ( Int_t j=0;j<TRACKR.mtrack;j++) {
b2be0e73 1488 sum +=TRACKR.dtrack[j];
829fb838 1489 }
1490 if (TRACKR.ntrack == 0 && TRACKR.mtrack == 0)
b2be0e73 1491 return fRull + sum;
829fb838 1492 else {
b2be0e73 1493 return sum;
829fb838 1494 }
1495}
1496
1497//______________________________________________________________________________
1498Int_t TFluka::TrackPid() const
1499{
1500// Return the id of the particle transported
1501// TRACKR.jtrack = identity number of the particle
a9ea1616 1502 FlukaCallerCode_t caller = GetCaller();
1503 if (caller != kEEDRAW) {
f926898e 1504 return PDGFromId(TRACKR.jtrack);
1505 }
829fb838 1506 else
1507 return -1000;
1508}
1509
1510//______________________________________________________________________________
1511Double_t TFluka::TrackCharge() const
1512{
1513// Return charge of the track currently transported
1514// PAPROP.ichrge = electric charge of the particle
1515// TRACKR.jtrack = identity number of the particle
a9ea1616 1516 FlukaCallerCode_t caller = GetCaller();
1517 if (caller != kEEDRAW)
829fb838 1518 return PAPROP.ichrge[TRACKR.jtrack+6];
1519 else
1520 return -1000.0;
1521}
1522
1523//______________________________________________________________________________
1524Double_t TFluka::TrackMass() const
1525{
1526// PAPROP.am = particle mass in GeV
1527// TRACKR.jtrack = identity number of the particle
a9ea1616 1528 FlukaCallerCode_t caller = GetCaller();
1529 if (caller != kEEDRAW)
829fb838 1530 return PAPROP.am[TRACKR.jtrack+6];
1531 else
1532 return -1000.0;
1533}
1534
1535//______________________________________________________________________________
1536Double_t TFluka::Etot() const
1537{
1538// TRACKR.etrack = total energy of the particle
a9ea1616 1539 FlukaCallerCode_t caller = GetCaller();
1540 if (caller != kEEDRAW)
829fb838 1541 return TRACKR.etrack;
1542 else
1543 return -1000.0;
1544}
1545
1546//
1547// track status
1548//
1549//______________________________________________________________________________
1550Bool_t TFluka::IsNewTrack() const
1551{
1552// Return true for the first call of Stepping()
1553 return fTrackIsNew;
1554}
1555
0dabe425 1556void TFluka::SetTrackIsNew(Bool_t flag)
1557{
1558// Return true for the first call of Stepping()
1559 fTrackIsNew = flag;
1560
1561}
1562
1563
829fb838 1564//______________________________________________________________________________
1565Bool_t TFluka::IsTrackInside() const
1566{
1567// True if the track is not at the boundary of the current volume
1568// In Fluka a step is always inside one kind of material
1569// If the step would go behind the region of one material,
1570// it will be shortened to reach only the boundary.
1571// Therefore IsTrackInside() is always true.
a9ea1616 1572 FlukaCallerCode_t caller = GetCaller();
1573 if (caller == kBXEntering || caller == kBXExiting)
829fb838 1574 return 0;
1575 else
1576 return 1;
1577}
1578
1579//______________________________________________________________________________
1580Bool_t TFluka::IsTrackEntering() const
1581{
1582// True if this is the first step of the track in the current volume
1583
a9ea1616 1584 FlukaCallerCode_t caller = GetCaller();
1585 if (caller == kBXEntering)
829fb838 1586 return 1;
1587 else return 0;
1588}
1589
1590//______________________________________________________________________________
1591Bool_t TFluka::IsTrackExiting() const
1592{
1593// True if track is exiting volume
1594//
a9ea1616 1595 FlukaCallerCode_t caller = GetCaller();
1596 if (caller == kBXExiting)
829fb838 1597 return 1;
1598 else return 0;
1599}
1600
1601//______________________________________________________________________________
1602Bool_t TFluka::IsTrackOut() const
1603{
1604// True if the track is out of the setup
1605// means escape
a9ea1616 1606 FlukaProcessCode_t icode = GetIcode();
1607
1608 if (icode == kKASKADescape ||
1609 icode == kEMFSCOescape ||
1610 icode == kKASNEUescape ||
1611 icode == kKASHEAescape ||
1612 icode == kKASOPHescape)
1613 return 1;
829fb838 1614 else return 0;
1615}
1616
1617//______________________________________________________________________________
1618Bool_t TFluka::IsTrackDisappeared() const
1619{
a9ea1616 1620// All inelastic interactions and decays
829fb838 1621// fIcode from usdraw
a9ea1616 1622 FlukaProcessCode_t icode = GetIcode();
1623 if (icode == kKASKADinelint || // inelastic interaction
1624 icode == kKASKADdecay || // particle decay
1625 icode == kKASKADdray || // delta ray generation by hadron
1626 icode == kKASKADpair || // direct pair production
1627 icode == kKASKADbrems || // bremsstrahlung (muon)
1628 icode == kEMFSCObrems || // bremsstrahlung (electron)
1629 icode == kEMFSCOmoller || // Moller scattering
1630 icode == kEMFSCObhabha || // Bhaba scattering
1631 icode == kEMFSCOanniflight || // in-flight annihilation
1632 icode == kEMFSCOannirest || // annihilation at rest
1633 icode == kEMFSCOpair || // pair production
1634 icode == kEMFSCOcompton || // Compton scattering
1635 icode == kEMFSCOphotoel || // Photoelectric effect
1636 icode == kKASNEUhadronic || // hadronic interaction
1637 icode == kKASHEAdray // delta-ray
0dabe425 1638 ) return 1;
829fb838 1639 else return 0;
1640}
1641
1642//______________________________________________________________________________
1643Bool_t TFluka::IsTrackStop() const
1644{
1645// True if the track energy has fallen below the threshold
1646// means stopped by signal or below energy threshold
a9ea1616 1647 FlukaProcessCode_t icode = GetIcode();
1648 if (icode == kKASKADstopping ||
1649 icode == kKASKADtimekill ||
1650 icode == kEMFSCOstopping1 ||
1651 icode == kEMFSCOstopping2 ||
1652 icode == kEMFSCOtimekill ||
1653 icode == kKASNEUstopping ||
1654 icode == kKASNEUtimekill ||
1655 icode == kKASHEAtimekill ||
1656 icode == kKASOPHtimekill) return 1;
829fb838 1657 else return 0;
1658}
1659
1660//______________________________________________________________________________
1661Bool_t TFluka::IsTrackAlive() const
1662{
1663// means not disappeared or not out
1664 if (IsTrackDisappeared() || IsTrackOut() ) return 0;
1665 else return 1;
1666}
1667
1668//
1669// secondaries
1670//
1671
1672//______________________________________________________________________________
1673Int_t TFluka::NSecondaries() const
1674
1675{
1676// Number of secondary particles generated in the current step
81f1d030 1677// GENSTK.np = number of secondaries except light and heavy ions
829fb838 1678// FHEAVY.npheav = number of secondaries for light and heavy secondary ions
a9ea1616 1679 FlukaCallerCode_t caller = GetCaller();
1680 if (caller == kUSDRAW) // valid only after usdraw
81f1d030 1681 return GENSTK.np + FHEAVY.npheav;
a9ea1616 1682 else if (caller == kUSTCKV) {
7b203b6e 1683 // Cerenkov Photon production
1684 return fNCerenkov;
1685 }
829fb838 1686 return 0;
1687} // end of NSecondaries
1688
1689//______________________________________________________________________________
1690void TFluka::GetSecondary(Int_t isec, Int_t& particleId,
1691 TLorentzVector& position, TLorentzVector& momentum)
1692{
1693// Copy particles from secondary stack to vmc stack
1694//
1695
a9ea1616 1696 FlukaCallerCode_t caller = GetCaller();
1697 if (caller == kUSDRAW) { // valid only after usdraw
81f1d030 1698 if (GENSTK.np > 0) {
7b203b6e 1699 // Hadronic interaction
81f1d030 1700 if (isec >= 0 && isec < GENSTK.np) {
1701 particleId = PDGFromId(GENSTK.kpart[isec]);
7b203b6e 1702 position.SetX(fXsco);
1703 position.SetY(fYsco);
1704 position.SetZ(fZsco);
1705 position.SetT(TRACKR.atrack);
81f1d030 1706 momentum.SetPx(GENSTK.plr[isec]*GENSTK.cxr[isec]);
1707 momentum.SetPy(GENSTK.plr[isec]*GENSTK.cyr[isec]);
1708 momentum.SetPz(GENSTK.plr[isec]*GENSTK.czr[isec]);
1709 momentum.SetE(GENSTK.tki[isec] + PAPROP.am[GENSTK.kpart[isec]+6]);
7b203b6e 1710 }
81f1d030 1711 else if (isec >= GENSTK.np && isec < GENSTK.np + FHEAVY.npheav) {
1712 Int_t jsec = isec - GENSTK.np;
7b203b6e 1713 particleId = FHEAVY.kheavy[jsec]; // this is Fluka id !!!
1714 position.SetX(fXsco);
1715 position.SetY(fYsco);
1716 position.SetZ(fZsco);
1717 position.SetT(TRACKR.atrack);
1718 momentum.SetPx(FHEAVY.pheavy[jsec]*FHEAVY.cxheav[jsec]);
1719 momentum.SetPy(FHEAVY.pheavy[jsec]*FHEAVY.cyheav[jsec]);
1720 momentum.SetPz(FHEAVY.pheavy[jsec]*FHEAVY.czheav[jsec]);
1721 if (FHEAVY.tkheav[jsec] >= 3 && FHEAVY.tkheav[jsec] <= 6)
1722 momentum.SetE(FHEAVY.tkheav[jsec] + PAPROP.am[jsec+6]);
1723 else if (FHEAVY.tkheav[jsec] > 6)
1724 momentum.SetE(FHEAVY.tkheav[jsec] + FHEAVY.amnhea[jsec]); // to be checked !!!
1725 }
1726 else
1727 Warning("GetSecondary","isec out of range");
1728 }
a9ea1616 1729 } else if (caller == kUSTCKV) {
7b203b6e 1730 Int_t index = OPPHST.lstopp - isec;
1731 position.SetX(OPPHST.xoptph[index]);
1732 position.SetY(OPPHST.yoptph[index]);
1733 position.SetZ(OPPHST.zoptph[index]);
1734 position.SetT(OPPHST.agopph[index]);
1735 Double_t p = OPPHST.poptph[index];
1736
1737 momentum.SetPx(p * OPPHST.txopph[index]);
1738 momentum.SetPy(p * OPPHST.tyopph[index]);
1739 momentum.SetPz(p * OPPHST.tzopph[index]);
1740 momentum.SetE(p);
829fb838 1741 }
1742 else
7b203b6e 1743 Warning("GetSecondary","no secondaries available");
1744
829fb838 1745} // end of GetSecondary
1746
7b203b6e 1747
829fb838 1748//______________________________________________________________________________
1749TMCProcess TFluka::ProdProcess(Int_t) const
1750
1751{
1752// Name of the process that has produced the secondary particles
1753// in the current step
0dabe425 1754
a9ea1616 1755 Int_t mugamma = (TRACKR.jtrack == kFLUKAphoton ||
1756 TRACKR.jtrack == kFLUKAmuplus ||
1757 TRACKR.jtrack == kFLUKAmuminus);
1758 FlukaProcessCode_t icode = GetIcode();
1759
1760 if (icode == kKASKADdecay) return kPDecay;
1761 else if (icode == kKASKADpair || icode == kEMFSCOpair) return kPPair;
1762 else if (icode == kEMFSCOcompton) return kPCompton;
1763 else if (icode == kEMFSCOphotoel) return kPPhotoelectric;
1764 else if (icode == kKASKADbrems || icode == kEMFSCObrems) return kPBrem;
1765 else if (icode == kKASKADdray || icode == kKASHEAdray) return kPDeltaRay;
1766 else if (icode == kEMFSCOmoller || icode == kEMFSCObhabha) return kPDeltaRay;
1767 else if (icode == kEMFSCOanniflight || icode == kEMFSCOannirest) return kPAnnihilation;
1768 else if (icode == kKASKADinelint) {
1769 if (!mugamma) return kPHadronic;
1770 else if (TRACKR.jtrack == kFLUKAphoton) return kPPhotoFission;
1771 else return kPMuonNuclear;
829fb838 1772 }
a9ea1616 1773 else if (icode == kEMFSCOrayleigh) return kPRayleigh;
829fb838 1774// Fluka codes 100, 300 and 400 still to be investigasted
a9ea1616 1775 else return kPNoProcess;
829fb838 1776}
1777
829fb838 1778
b496f27c 1779Int_t TFluka::StepProcesses(TArrayI &proc) const
1780{
1781 //
1782 // Return processes active in the current step
1783 //
a9ea1616 1784 FlukaProcessCode_t icode = GetIcode();
b496f27c 1785 proc.Set(1);
1786 TMCProcess iproc;
a9ea1616 1787 switch (icode) {
1788 case kKASKADtimekill:
1789 case kEMFSCOtimekill:
1790 case kKASNEUtimekill:
1791 case kKASHEAtimekill:
1792 case kKASOPHtimekill:
b496f27c 1793 iproc = kPTOFlimit;
1794 break;
a9ea1616 1795 case kKASKADstopping:
1796 case kKASKADescape:
1797 case kEMFSCOstopping1:
1798 case kEMFSCOstopping2:
1799 case kEMFSCOescape:
1800 case kKASNEUstopping:
1801 case kKASNEUescape:
1802 case kKASHEAescape:
1803 case kKASOPHescape:
6fd5baa4 1804 iproc = kPStop;
b496f27c 1805 break;
a9ea1616 1806 case kKASOPHabsorption:
b496f27c 1807 iproc = kPLightAbsorption;
1808 break;
a9ea1616 1809 case kKASOPHrefraction:
6fd5baa4 1810 iproc = kPLightRefraction;
a9ea1616 1811 case kEMSCOlocaledep :
b496f27c 1812 iproc = kPPhotoelectric;
1813 break;
1814 default:
1815 iproc = ProdProcess(0);
1816 }
07f5b33e 1817 proc[0] = iproc;
b496f27c 1818 return 1;
1819}
829fb838 1820//______________________________________________________________________________
1821Int_t TFluka::VolId2Mate(Int_t id) const
1822{
1823//
1824// Returns the material number for a given volume ID
1825//
1826 return fMCGeo->VolId2Mate(id);
1827}
1828
1829//______________________________________________________________________________
1830const char* TFluka::VolName(Int_t id) const
1831{
1832//
1833// Returns the volume name for a given volume ID
1834//
1835 return fMCGeo->VolName(id);
1836}
1837
1838//______________________________________________________________________________
1839Int_t TFluka::VolId(const Text_t* volName) const
1840{
1841//
1842// Converts from volume name to volume ID.
1843// Time consuming. (Only used during set-up)
1844// Could be replaced by hash-table
1845//
09cd6497 1846 char sname[20];
1847 Int_t len;
1848 strncpy(sname, volName, len = strlen(volName));
1849 sname[len] = 0;
1850 while (sname[len - 1] == ' ') sname[--len] = 0;
1851 return fMCGeo->VolId(sname);
829fb838 1852}
1853
1854//______________________________________________________________________________
1855Int_t TFluka::CurrentVolID(Int_t& copyNo) const
1856{
1857//
1858// Return the logical id and copy number corresponding to the current fluka region
1859//
1860 if (gGeoManager->IsOutside()) return 0;
1861 TGeoNode *node = gGeoManager->GetCurrentNode();
1862 copyNo = node->GetNumber();
1863 Int_t id = node->GetVolume()->GetNumber();
1864 return id;
1865}
1866
1867//______________________________________________________________________________
1868Int_t TFluka::CurrentVolOffID(Int_t off, Int_t& copyNo) const
1869{
1870//
1871// Return the logical id and copy number of off'th mother
1872// corresponding to the current fluka region
1873//
1874 if (off<0 || off>gGeoManager->GetLevel()) return 0;
1875 if (off==0) return CurrentVolID(copyNo);
1876 TGeoNode *node = gGeoManager->GetMother(off);
1877 if (!node) return 0;
1878 copyNo = node->GetNumber();
1879 return node->GetVolume()->GetNumber();
1880}
1881
1882//______________________________________________________________________________
1883const char* TFluka::CurrentVolName() const
1884{
1885//
1886// Return the current volume name
1887//
1888 if (gGeoManager->IsOutside()) return 0;
1889 return gGeoManager->GetCurrentVolume()->GetName();
1890}
1891
1892//______________________________________________________________________________
1893const char* TFluka::CurrentVolOffName(Int_t off) const
1894{
1895//
1896// Return the volume name of the off'th mother of the current volume
1897//
1898 if (off<0 || off>gGeoManager->GetLevel()) return 0;
1899 if (off==0) return CurrentVolName();
1900 TGeoNode *node = gGeoManager->GetMother(off);
1901 if (!node) return 0;
1902 return node->GetVolume()->GetName();
1903}
1904
d59acfe7 1905const char* TFluka::CurrentVolPath() {
1906 // Return the current volume path
1907 return gGeoManager->GetPath();
1908}
829fb838 1909//______________________________________________________________________________
a60813de 1910Int_t TFluka::CurrentMaterial(Float_t & a, Float_t & z,
1911 Float_t & dens, Float_t & radl, Float_t & absl) const
829fb838 1912{
1913//
a60813de 1914// Return the current medium number and material properties
829fb838 1915//
1916 Int_t copy;
1917 Int_t id = TFluka::CurrentVolID(copy);
1918 Int_t med = TFluka::VolId2Mate(id);
a60813de 1919 TGeoVolume* vol = gGeoManager->GetCurrentVolume();
1920 TGeoMaterial* mat = vol->GetMaterial();
1921 a = mat->GetA();
1922 z = mat->GetZ();
1923 dens = mat->GetDensity();
1924 radl = mat->GetRadLen();
1925 absl = mat->GetIntLen();
1926
829fb838 1927 return med;
1928}
1929
1930//______________________________________________________________________________
1931void TFluka::Gmtod(Float_t* xm, Float_t* xd, Int_t iflag)
1932{
1933// Transforms a position from the world reference frame
1934// to the current volume reference frame.
1935//
1936// Geant3 desription:
1937// ==================
1938// Computes coordinates XD (in DRS)
1939// from known coordinates XM in MRS
1940// The local reference system can be initialized by
1941// - the tracking routines and GMTOD used in GUSTEP
1942// - a call to GMEDIA(XM,NUMED)
1943// - a call to GLVOLU(NLEVEL,NAMES,NUMBER,IER)
1944// (inverse routine is GDTOM)
1945//
1946// If IFLAG=1 convert coordinates
1947// IFLAG=2 convert direction cosinus
1948//
1949// ---
1950 Double_t xmL[3], xdL[3];
1951 Int_t i;
1952 for (i=0;i<3;i++) xmL[i]=xm[i];
1953 if (iflag == 1) gGeoManager->MasterToLocal(xmL,xdL);
1954 else gGeoManager->MasterToLocalVect(xmL,xdL);
1955 for (i=0;i<3;i++) xd[i] = xdL[i];
1956}
1957
1958//______________________________________________________________________________
1959void TFluka::Gmtod(Double_t* xm, Double_t* xd, Int_t iflag)
1960{
1961 if (iflag == 1) gGeoManager->MasterToLocal(xm,xd);
1962 else gGeoManager->MasterToLocalVect(xm,xd);
1963}
1964
1965//______________________________________________________________________________
1966void TFluka::Gdtom(Float_t* xd, Float_t* xm, Int_t iflag)
1967{
1968// Transforms a position from the current volume reference frame
1969// to the world reference frame.
1970//
1971// Geant3 desription:
1972// ==================
1973// Computes coordinates XM (Master Reference System
1974// knowing the coordinates XD (Detector Ref System)
1975// The local reference system can be initialized by
1976// - the tracking routines and GDTOM used in GUSTEP
1977// - a call to GSCMED(NLEVEL,NAMES,NUMBER)
1978// (inverse routine is GMTOD)
1979//
1980// If IFLAG=1 convert coordinates
1981// IFLAG=2 convert direction cosinus
1982//
1983// ---
1984 Double_t xmL[3], xdL[3];
1985 Int_t i;
1986 for (i=0;i<3;i++) xdL[i] = xd[i];
1987 if (iflag == 1) gGeoManager->LocalToMaster(xdL,xmL);
1988 else gGeoManager->LocalToMasterVect(xdL,xmL);
1989 for (i=0;i<3;i++) xm[i]=xmL[i];
1990}
1991
1992//______________________________________________________________________________
1993void TFluka::Gdtom(Double_t* xd, Double_t* xm, Int_t iflag)
1994{
1995 if (iflag == 1) gGeoManager->LocalToMaster(xd,xm);
1996 else gGeoManager->LocalToMasterVect(xd,xm);
1997}
1998
1999//______________________________________________________________________________
2000TObjArray *TFluka::GetFlukaMaterials()
2001{
2002 return fGeom->GetMatList();
2003}
2004
2005//______________________________________________________________________________
a9ea1616 2006void TFluka::SetMreg(Int_t l, Int_t lttc)
829fb838 2007{
2008// Set current fluka region
2009 fCurrentFlukaRegion = l;
a9ea1616 2010 fGeom->SetMreg(l,lttc);
829fb838 2011}
2012
2013
b496f27c 2014
2015
2016TString TFluka::ParticleName(Int_t pdg) const
2017{
2018 // Return particle name for particle with pdg code pdg.
2019 Int_t ifluka = IdFromPDG(pdg);
ece92b30 2020 return TString((CHPPRP.btype[ifluka - kFLUKAcodemin]), 8);
b496f27c 2021}
2022
2023
2024Double_t TFluka::ParticleMass(Int_t pdg) const
2025{
2026 // Return particle mass for particle with pdg code pdg.
2027 Int_t ifluka = IdFromPDG(pdg);
ece92b30 2028 return (PAPROP.am[ifluka - kFLUKAcodemin]);
2029}
2030
2031Double_t TFluka::ParticleMassFPC(Int_t fpc) const
2032{
2033 // Return particle mass for particle with Fluka particle code fpc
2034 return (PAPROP.am[fpc - kFLUKAcodemin]);
b496f27c 2035}
2036
2037Double_t TFluka::ParticleCharge(Int_t pdg) const
2038{
2039 // Return particle charge for particle with pdg code pdg.
2040 Int_t ifluka = IdFromPDG(pdg);
ece92b30 2041 return Double_t(PAPROP.ichrge[ifluka - kFLUKAcodemin]);
b496f27c 2042}
2043
2044Double_t TFluka::ParticleLifeTime(Int_t pdg) const
2045{
2046 // Return particle lifetime for particle with pdg code pdg.
2047 Int_t ifluka = IdFromPDG(pdg);
ece92b30 2048 return (PAPROP.tmnlf[ifluka - kFLUKAcodemin]);
b496f27c 2049}
2050
2051void TFluka::Gfpart(Int_t pdg, char* name, Int_t& type, Float_t& mass, Float_t& charge, Float_t& tlife)
2052{
2053 // Retrieve particle properties for particle with pdg code pdg.
2054
2055 strcpy(name, ParticleName(pdg).Data());
2056 type = ParticleMCType(pdg);
2057 mass = ParticleMass(pdg);
2058 charge = ParticleCharge(pdg);
2059 tlife = ParticleLifeTime(pdg);
2060}
2061
8e5bf079 2062void TFluka::PrintHeader()
2063{
2064 //
2065 // Print a header
2066 printf("\n");
2067 printf("\n");
2068 printf("------------------------------------------------------------------------------\n");
2069 printf("- You are using the TFluka Virtual Monte Carlo Interface to FLUKA. -\n");
2070 printf("- Please see the file fluka.out for FLUKA output and licensing information. -\n");
2071 printf("------------------------------------------------------------------------------\n");
2072 printf("\n");
2073 printf("\n");
2074}
2075
b496f27c 2076
81f1d030 2077#define pshckp pshckp_
2078#define ustckv ustckv_
3a625972 2079
2080
2081extern "C" {
81f1d030 2082 void pshckp(Double_t & px, Double_t & py, Double_t & pz, Double_t & e,
2083 Double_t & vx, Double_t & vy, Double_t & vz, Double_t & tof,
2084 Double_t & polx, Double_t & poly, Double_t & polz, Double_t & wgt, Int_t& ntr)
2085 {
2086 //
2087 // Pushes one cerenkov photon to the stack
2088 //
2089
2090 TFluka* fluka = (TFluka*) gMC;
2091 TVirtualMCStack* cppstack = fluka->GetStack();
2092 Int_t parent = TRACKR.ispusr[mkbmx2-1];
2093 cppstack->PushTrack(0, parent, 50000050,
2094 px, py, pz, e,
2095 vx, vy, vz, tof,
2096 polx, poly, polz,
2097 kPCerenkov, ntr, wgt, 0);
2098 }
2099
2100 void ustckv(Int_t & nphot, Int_t & mreg, Double_t & x, Double_t & y, Double_t & z)
7b203b6e 2101 {
2102 //
2103 // Calls stepping in order to signal cerenkov production
2104 //
2105 TFluka *fluka = (TFluka*)gMC;
a9ea1616 2106 fluka->SetMreg(mreg,LTCLCM.mlatm1);
7b203b6e 2107 fluka->SetXsco(x);
2108 fluka->SetYsco(y);
2109 fluka->SetZsco(z);
2110 fluka->SetNCerenkov(nphot);
a9ea1616 2111 fluka->SetCaller(kUSTCKV);
d59acfe7 2112 if (fluka->GetVerbosityLevel() >= 3)
7b203b6e 2113 (TVirtualMCApplication::Instance())->Stepping();
81f1d030 2114
7b203b6e 2115 }
3a625972 2116}
a9ea1616 2117