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