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