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