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