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