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