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