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