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