Removing objects in a TClonesArray
[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"
35#include "TCallf77.h" //For the fortran calls
36#include "Fdblprc.h" //(DBLPRC) fluka common
37#include "Fepisor.h" //(EPISOR) fluka common
07f5b33e 38#include "Ffinuc.h" //(FINUC) fluka common
829fb838 39#include "Fiounit.h" //(IOUNIT) fluka common
40#include "Fpaprop.h" //(PAPROP) fluka common
41#include "Fpart.h" //(PART) fluka common
42#include "Ftrackr.h" //(TRACKR) fluka common
43#include "Fpaprop.h" //(PAPROP) fluka common
44#include "Ffheavy.h" //(FHEAVY) fluka common
3a625972 45#include "Fopphst.h" //(OPPHST) fluka common
07f5b33e 46#include "Fstack.h" //(STACK) fluka common
47#include "Fstepsz.h" //(STEPSZ) fluka common
7b203b6e 48#include "Fopphst.h" //(OPPHST) fluka common
829fb838 49
50#include "TVirtualMC.h"
3a625972 51#include "TMCProcess.h"
829fb838 52#include "TGeoManager.h"
53#include "TGeoMaterial.h"
54#include "TGeoMedium.h"
55#include "TFlukaMCGeometry.h"
6f5667d1 56#include "TGeoMCGeometry.h"
829fb838 57#include "TFlukaCerenkov.h"
1df5fa54 58#include "TFlukaConfigOption.h"
b496f27c 59#include "TFlukaScoringOption.h"
829fb838 60#include "TLorentzVector.h"
b496f27c 61#include "TArrayI.h"
829fb838 62
63// Fluka methods that may be needed.
64#ifndef WIN32
65# define flukam flukam_
66# define fluka_openinp fluka_openinp_
67# define fluka_closeinp fluka_closeinp_
68# define mcihad mcihad_
69# define mpdgha mpdgha_
eea53470 70# define newplo newplo_
829fb838 71#else
72# define flukam FLUKAM
73# define fluka_openinp FLUKA_OPENINP
74# define fluka_closeinp FLUKA_CLOSEINP
75# define mcihad MCIHAD
76# define mpdgha MPDGHA
eea53470 77# define newplo NEWPLO
829fb838 78#endif
79
80extern "C"
81{
82 //
83 // Prototypes for FLUKA functions
84 //
85 void type_of_call flukam(const int&);
eea53470 86 void type_of_call newplo();
829fb838 87 void type_of_call fluka_openinp(const int&, DEFCHARA);
88 void type_of_call fluka_closeinp(const int&);
89 int type_of_call mcihad(const int&);
90 int type_of_call mpdgha(const int&);
91}
92
93//
94// Class implementation for ROOT
95//
96ClassImp(TFluka)
97
98//
99//----------------------------------------------------------------------------
100// TFluka constructors and destructors.
101//______________________________________________________________________________
102TFluka::TFluka()
103 :TVirtualMC(),
104 fVerbosityLevel(0),
1df5fa54 105 fInputFileName(""),
fb2cbbec 106 fUserConfig(0),
1df5fa54 107 fUserScore(0)
829fb838 108{
109 //
110 // Default constructor
111 //
112 fGeneratePemf = kFALSE;
113 fNVolumes = 0;
114 fCurrentFlukaRegion = -1;
115 fGeom = 0;
116 fMCGeo = 0;
117 fMaterials = 0;
118 fDummyBoundary = 0;
119 fFieldFlag = 1;
bd3d5c8a 120 fStopped = 0;
b496f27c 121 fStopEvent = 0;
122 fStopRun = 0;
123 fNEvent = 0;
829fb838 124}
125
126//______________________________________________________________________________
127TFluka::TFluka(const char *title, Int_t verbosity, Bool_t isRootGeometrySupported)
128 :TVirtualMC("TFluka",title, isRootGeometrySupported),
129 fVerbosityLevel(verbosity),
130 fInputFileName(""),
131 fTrackIsEntering(0),
132 fTrackIsExiting(0),
1df5fa54 133 fTrackIsNew(0),
fb2cbbec 134 fUserConfig(new TObjArray(100)),
1df5fa54 135 fUserScore(new TObjArray(100))
829fb838 136{
137 // create geometry interface
7f13be31 138 if (fVerbosityLevel >=3)
139 cout << "<== TFluka::TFluka(" << title << ") constructor called." << endl;
140 SetCoreInputFileName();
141 SetInputFileName();
142 SetGeneratePemf(kFALSE);
829fb838 143 fNVolumes = 0;
144 fCurrentFlukaRegion = -1;
145 fDummyBoundary = 0;
146 fFieldFlag = 1;
147 fGeneratePemf = kFALSE;
148 fMCGeo = new TGeoMCGeometry("MCGeo", "TGeo Implementation of VirtualMCGeometry", kTRUE);
fb2cbbec 149 fGeom = new TFlukaMCGeometry("geom", "FLUKA VMC Geometry");
829fb838 150 if (verbosity > 2) fGeom->SetDebugMode(kTRUE);
151 fMaterials = 0;
bd3d5c8a 152 fStopped = 0;
b496f27c 153 fStopEvent = 0;
154 fStopRun = 0;
155 fNEvent = 0;
829fb838 156}
157
158//______________________________________________________________________________
159TFluka::~TFluka() {
160// Destructor
1df5fa54 161 if (fVerbosityLevel >=3)
162 cout << "<== TFluka::~TFluka() destructor called." << endl;
163
164 delete fGeom;
165 delete fMCGeo;
166
fb2cbbec 167 if (fUserConfig) {
168 fUserConfig->Delete();
169 delete fUserConfig;
1df5fa54 170 }
171
172
829fb838 173}
174
175//
176//______________________________________________________________________________
177// TFluka control methods
178//______________________________________________________________________________
179void TFluka::Init() {
180//
181// Geometry initialisation
182//
183 if (fVerbosityLevel >=3) cout << "==> TFluka::Init() called." << endl;
184
185 if (!gGeoManager) new TGeoManager("geom", "FLUKA geometry");
186 fApplication->ConstructGeometry();
187 TGeoVolume *top = (TGeoVolume*)gGeoManager->GetListOfVolumes()->First();
188 gGeoManager->SetTopVolume(top);
189 gGeoManager->CloseGeometry("di");
190 gGeoManager->DefaultColors(); // to be removed
191 fNVolumes = fGeom->NofVolumes();
192 fGeom->CreateFlukaMatFile("flukaMat.inp");
193 if (fVerbosityLevel >=3) {
194 printf("== Number of volumes: %i\n ==", fNVolumes);
195 cout << "\t* InitPhysics() - Prepare input file to be called" << endl;
196 }
3b8c325d 197 // now we have TGeo geometry created and we have to patch FlukaVmc.inp
829fb838 198 // with the material mapping file FlukaMat.inp
199}
200
201
202//______________________________________________________________________________
203void TFluka::FinishGeometry() {
204//
205// Build-up table with region to medium correspondance
206//
207 if (fVerbosityLevel >=3) {
208 cout << "==> TFluka::FinishGeometry() called." << endl;
209 printf("----FinishGeometry - nothing to do with TGeo\n");
210 cout << "<== TFluka::FinishGeometry() called." << endl;
211 }
212}
213
214//______________________________________________________________________________
215void TFluka::BuildPhysics() {
216//
217// Prepare FLUKA input files and call FLUKA physics initialisation
218//
219
220 if (fVerbosityLevel >=3)
221 cout << "==> TFluka::BuildPhysics() called." << endl;
222// Prepare input file with the current physics settings
223 InitPhysics();
224 cout << "\t* InitPhysics() - Prepare input file was called" << endl;
225
226 if (fVerbosityLevel >=2)
227 cout << "\t* Changing lfdrtr = (" << (GLOBAL.lfdrtr?'T':'F')
228 << ") in fluka..." << endl;
229 GLOBAL.lfdrtr = true;
230
231 if (fVerbosityLevel >=2)
232 cout << "\t* Opening file " << fInputFileName << endl;
233 const char* fname = fInputFileName;
234 fluka_openinp(lunin, PASSCHARA(fname));
235
236 if (fVerbosityLevel >=2)
237 cout << "\t* Calling flukam..." << endl;
238 flukam(1);
239
240 if (fVerbosityLevel >=2)
241 cout << "\t* Closing file " << fInputFileName << endl;
242 fluka_closeinp(lunin);
243
244 FinishGeometry();
245
246 if (fVerbosityLevel >=3)
247 cout << "<== TFluka::Init() called." << endl;
248
249
250 if (fVerbosityLevel >=3)
251 cout << "<== TFluka::BuildPhysics() called." << endl;
252}
253
254//______________________________________________________________________________
255void TFluka::ProcessEvent() {
256//
257// Process one event
258//
b496f27c 259 if (fStopRun) {
260 printf("User Run Abortion: No more events handled !\n");
261 fNEvent += 1;
262 return;
263 }
264
265 if (fVerbosityLevel >=3)
266 cout << "==> TFluka::ProcessEvent() called." << endl;
267 fApplication->GeneratePrimaries();
268 EPISOR.lsouit = true;
269 flukam(1);
270 if (fVerbosityLevel >=3)
271 cout << "<== TFluka::ProcessEvent() called." << endl;
272 //
273 // Increase event number
274 //
275 fNEvent += 1;
829fb838 276}
277
278//______________________________________________________________________________
279Bool_t TFluka::ProcessRun(Int_t nevent) {
280//
281// Run steering
282//
283
284 if (fVerbosityLevel >=3)
285 cout << "==> TFluka::ProcessRun(" << nevent << ") called."
286 << endl;
287
288 if (fVerbosityLevel >=2) {
289 cout << "\t* GLOBAL.fdrtr = " << (GLOBAL.lfdrtr?'T':'F') << endl;
290 cout << "\t* Calling flukam again..." << endl;
291 }
292
293 fApplication->InitGeometry();
294 Int_t todo = TMath::Abs(nevent);
295 for (Int_t ev = 0; ev < todo; ev++) {
296 fApplication->BeginEvent();
297 ProcessEvent();
298 fApplication->FinishEvent();
299 }
300
301 if (fVerbosityLevel >=3)
302 cout << "<== TFluka::ProcessRun(" << nevent << ") called."
303 << endl;
eea53470 304 // Write fluka specific scoring output
305 newplo();
306
829fb838 307 return kTRUE;
308}
309
310//_____________________________________________________________________________
311// methods for building/management of geometry
312
313// functions from GCONS
314//____________________________________________________________________________
315void TFluka::Gfmate(Int_t imat, char *name, Float_t &a, Float_t &z,
316 Float_t &dens, Float_t &radl, Float_t &absl,
317 Float_t* /*ubuf*/, Int_t& /*nbuf*/) {
318//
319 TGeoMaterial *mat;
320 TIter next (gGeoManager->GetListOfMaterials());
321 while ((mat = (TGeoMaterial*)next())) {
322 if (mat->GetUniqueID() == (UInt_t)imat) break;
323 }
324 if (!mat) {
325 Error("Gfmate", "no material with index %i found", imat);
326 return;
327 }
328 sprintf(name, "%s", mat->GetName());
329 a = mat->GetA();
330 z = mat->GetZ();
331 dens = mat->GetDensity();
332 radl = mat->GetRadLen();
333 absl = mat->GetIntLen();
334}
335
336//______________________________________________________________________________
337void TFluka::Gfmate(Int_t imat, char *name, Double_t &a, Double_t &z,
338 Double_t &dens, Double_t &radl, Double_t &absl,
339 Double_t* /*ubuf*/, Int_t& /*nbuf*/) {
340//
341 TGeoMaterial *mat;
342 TIter next (gGeoManager->GetListOfMaterials());
343 while ((mat = (TGeoMaterial*)next())) {
344 if (mat->GetUniqueID() == (UInt_t)imat) break;
345 }
346 if (!mat) {
347 Error("Gfmate", "no material with index %i found", imat);
348 return;
349 }
350 sprintf(name, "%s", mat->GetName());
351 a = mat->GetA();
352 z = mat->GetZ();
353 dens = mat->GetDensity();
354 radl = mat->GetRadLen();
355 absl = mat->GetIntLen();
356}
357
358// detector composition
359//______________________________________________________________________________
360void TFluka::Material(Int_t& kmat, const char* name, Double_t a,
361 Double_t z, Double_t dens, Double_t radl, Double_t absl,
362 Float_t* buf, Int_t nwbuf) {
363//
364 Double_t* dbuf = fGeom->CreateDoubleArray(buf, nwbuf);
365 Material(kmat, name, a, z, dens, radl, absl, dbuf, nwbuf);
366 delete [] dbuf;
367}
368
369//______________________________________________________________________________
370void TFluka::Material(Int_t& kmat, const char* name, Double_t a,
371 Double_t z, Double_t dens, Double_t radl, Double_t absl,
372 Double_t* /*buf*/, Int_t /*nwbuf*/) {
373//
fb2cbbec 374// Define a material
829fb838 375 TGeoMaterial *mat;
376 kmat = gGeoManager->GetListOfMaterials()->GetSize();
377 if ((z-Int_t(z)) > 1E-3) {
378 mat = fGeom->GetMakeWrongMaterial(z);
379 if (mat) {
380 mat->SetRadLen(radl,absl);
381 mat->SetUniqueID(kmat);
382 return;
383 }
384 }
385 gGeoManager->Material(name, a, z, dens, kmat, radl, absl);
386}
387
388//______________________________________________________________________________
389void TFluka::Mixture(Int_t& kmat, const char *name, Float_t *a,
390 Float_t *z, Double_t dens, Int_t nlmat, Float_t *wmat) {
391//
fb2cbbec 392// Define a material mixture
393//
829fb838 394 Double_t* da = fGeom->CreateDoubleArray(a, TMath::Abs(nlmat));
395 Double_t* dz = fGeom->CreateDoubleArray(z, TMath::Abs(nlmat));
396 Double_t* dwmat = fGeom->CreateDoubleArray(wmat, TMath::Abs(nlmat));
397
398 Mixture(kmat, name, da, dz, dens, nlmat, dwmat);
399 for (Int_t i=0; i<nlmat; i++) {
400 a[i] = da[i]; z[i] = dz[i]; wmat[i] = dwmat[i];
401 }
402
403 delete [] da;
404 delete [] dz;
405 delete [] dwmat;
406}
407
408//______________________________________________________________________________
409void TFluka::Mixture(Int_t& kmat, const char *name, Double_t *a,
410 Double_t *z, Double_t dens, Int_t nlmat, Double_t *wmat) {
411//
412 // Defines mixture OR COMPOUND IMAT as composed by
413 // THE BASIC NLMAT materials defined by arrays A,Z and WMAT
414 //
415 // If NLMAT > 0 then wmat contains the proportion by
416 // weights of each basic material in the mixture.
417 //
418 // If nlmat < 0 then WMAT contains the number of atoms
419 // of a given kind into the molecule of the COMPOUND
420 // In this case, WMAT in output is changed to relative
421 // weigths.
422 //
423 Int_t i,j;
424 if (nlmat < 0) {
425 nlmat = - nlmat;
426 Double_t amol = 0;
427 for (i=0;i<nlmat;i++) {
428 amol += a[i]*wmat[i];
429 }
430 for (i=0;i<nlmat;i++) {
431 wmat[i] *= a[i]/amol;
432 }
433 }
434 kmat = gGeoManager->GetListOfMaterials()->GetSize();
435 // Check if we have elements with fractional Z
436 TGeoMaterial *mat = 0;
437 TGeoMixture *mix = 0;
438 Bool_t mixnew = kFALSE;
439 for (i=0; i<nlmat; i++) {
440 if (z[i]-Int_t(z[i]) < 1E-3) continue;
441 // We have found an element with fractional Z -> loop mixtures to look for it
442 for (j=0; j<kmat; j++) {
443 mat = (TGeoMaterial*)gGeoManager->GetListOfMaterials()->At(j);
444 if (!mat) break;
445 if (!mat->IsMixture()) continue;
446 mix = (TGeoMixture*)mat;
447 if (TMath::Abs(z[i]-mix->GetZ()) >1E-3) continue;
448// printf(" FOUND component %i as mixture %s\n", i, mat->GetName());
449 mixnew = kTRUE;
450 break;
451 }
452 if (!mixnew) Warning("Mixture","%s : cannot find component %i with fractional Z=%f\n", name, i, z[i]);
453 break;
454 }
455 if (mixnew) {
456 Int_t nlmatnew = nlmat+mix->GetNelements()-1;
457 Double_t *anew = new Double_t[nlmatnew];
458 Double_t *znew = new Double_t[nlmatnew];
459 Double_t *wmatnew = new Double_t[nlmatnew];
460 Int_t ind=0;
461 for (j=0; j<nlmat; j++) {
462 if (j==i) continue;
463 anew[ind] = a[j];
464 znew[ind] = z[j];
465 wmatnew[ind] = wmat[j];
466 ind++;
467 }
468 for (j=0; j<mix->GetNelements(); j++) {
469 anew[ind] = mix->GetAmixt()[j];
470 znew[ind] = mix->GetZmixt()[j];
471 wmatnew[ind] = wmat[i]*mix->GetWmixt()[j];
472 ind++;
473 }
474 Mixture(kmat, name, anew, znew, dens, nlmatnew, wmatnew);
475 delete [] anew;
476 delete [] znew;
477 delete [] wmatnew;
478 return;
479 }
480 // Now we need to compact identical elements within the mixture
481 // First check if this happens
482 mixnew = kFALSE;
483 for (i=0; i<nlmat-1; i++) {
484 for (j=i+1; j<nlmat; j++) {
485 if (z[i] == z[j]) {
486 mixnew = kTRUE;
487 break;
488 }
489 }
490 if (mixnew) break;
491 }
492 if (mixnew) {
493 Int_t nlmatnew = 0;
494 Double_t *anew = new Double_t[nlmat];
495 Double_t *znew = new Double_t[nlmat];
496 memset(znew, 0, nlmat*sizeof(Double_t));
497 Double_t *wmatnew = new Double_t[nlmat];
498 Bool_t skipi;
499 for (i=0; i<nlmat; i++) {
500 skipi = kFALSE;
501 for (j=0; j<nlmatnew; j++) {
502 if (z[i] == z[j]) {
503 wmatnew[j] += wmat[i];
504 skipi = kTRUE;
505 break;
506 }
507 }
508 if (skipi) continue;
509 anew[nlmatnew] = a[i];
510 znew[nlmatnew] = z[i];
511 wmatnew[nlmatnew] = wmat[i];
512 nlmatnew++;
513 }
514 Mixture(kmat, name, anew, znew, dens, nlmatnew, wmatnew);
515 delete [] anew;
516 delete [] znew;
517 delete [] wmatnew;
518 return;
519 }
520 gGeoManager->Mixture(name, a, z, dens, nlmat, wmat, kmat);
521}
522
523//______________________________________________________________________________
524void TFluka::Medium(Int_t& kmed, const char *name, Int_t nmat,
525 Int_t isvol, Int_t ifield, Double_t fieldm, Double_t tmaxfd,
526 Double_t stemax, Double_t deemax, Double_t epsil,
527 Double_t stmin, Float_t* ubuf, Int_t nbuf) {
b2129742 528 // Define a medium
529 //
829fb838 530 kmed = gGeoManager->GetListOfMedia()->GetSize()+1;
531 fMCGeo->Medium(kmed, name, nmat, isvol, ifield, fieldm, tmaxfd, stemax, deemax,
532 epsil, stmin, ubuf, nbuf);
533}
534
535//______________________________________________________________________________
536void TFluka::Medium(Int_t& kmed, const char *name, Int_t nmat,
537 Int_t isvol, Int_t ifield, Double_t fieldm, Double_t tmaxfd,
538 Double_t stemax, Double_t deemax, Double_t epsil,
539 Double_t stmin, Double_t* ubuf, Int_t nbuf) {
b2129742 540 // Define a medium
541 //
829fb838 542 kmed = gGeoManager->GetListOfMedia()->GetSize()+1;
543 fMCGeo->Medium(kmed, name, nmat, isvol, ifield, fieldm, tmaxfd, stemax, deemax,
544 epsil, stmin, ubuf, nbuf);
545}
546
547//______________________________________________________________________________
548void TFluka::Matrix(Int_t& krot, Double_t thetaX, Double_t phiX,
549 Double_t thetaY, Double_t phiY, Double_t thetaZ,
550 Double_t phiZ) {
551//
552 krot = gGeoManager->GetListOfMatrices()->GetEntriesFast();
553 fMCGeo->Matrix(krot, thetaX, phiX, thetaY, phiY, thetaZ, phiZ);
554}
555
556//______________________________________________________________________________
557void TFluka::Gstpar(Int_t itmed, const char* param, Double_t parval) {
558//
559//
7b203b6e 560// Check if material is used
561 if (fVerbosityLevel >=3)
562 printf("Gstpar called with %6d %5s %12.4e %6d\n", itmed, param, parval, fGeom->GetFlukaMaterial(itmed));
563 Int_t* reglist;
564 Int_t nreg;
c1c801f9 565 reglist = fGeom->GetMaterialList(fGeom->GetFlukaMaterial(itmed), nreg);
566 if (nreg == 0) {
567 return;
568 }
569
7b203b6e 570//
829fb838 571 Bool_t process = kFALSE;
572 if (strncmp(param, "DCAY", 4) == 0 ||
573 strncmp(param, "PAIR", 4) == 0 ||
574 strncmp(param, "COMP", 4) == 0 ||
575 strncmp(param, "PHOT", 4) == 0 ||
576 strncmp(param, "PFIS", 4) == 0 ||
577 strncmp(param, "DRAY", 4) == 0 ||
578 strncmp(param, "ANNI", 4) == 0 ||
579 strncmp(param, "BREM", 4) == 0 ||
580 strncmp(param, "MUNU", 4) == 0 ||
581 strncmp(param, "CKOV", 4) == 0 ||
582 strncmp(param, "HADR", 4) == 0 ||
583 strncmp(param, "LOSS", 4) == 0 ||
584 strncmp(param, "MULS", 4) == 0 ||
585 strncmp(param, "RAYL", 4) == 0)
586 {
587 process = kTRUE;
588 }
589 if (process) {
590 SetProcess(param, Int_t (parval), fGeom->GetFlukaMaterial(itmed));
591 } else {
592 SetCut(param, parval, fGeom->GetFlukaMaterial(itmed));
593 }
594}
595
596// functions from GGEOM
597//_____________________________________________________________________________
598void TFluka::Gsatt(const char *name, const char *att, Int_t val)
599{
6f5667d1 600 // Set visualisation attributes for one volume
829fb838 601 char vname[5];
602 fGeom->Vname(name,vname);
603 char vatt[5];
604 fGeom->Vname(att,vatt);
605 gGeoManager->SetVolumeAttribute(vname, vatt, val);
606}
607
608//______________________________________________________________________________
609Int_t TFluka::Gsvolu(const char *name, const char *shape, Int_t nmed,
610 Float_t *upar, Int_t np) {
611//
612 return fMCGeo->Gsvolu(name, shape, nmed, upar, np);
613}
614
615//______________________________________________________________________________
616Int_t TFluka::Gsvolu(const char *name, const char *shape, Int_t nmed,
617 Double_t *upar, Int_t np) {
618//
619 return fMCGeo->Gsvolu(name, shape, nmed, upar, np);
620}
621
622//______________________________________________________________________________
623void TFluka::Gsdvn(const char *name, const char *mother, Int_t ndiv,
624 Int_t iaxis) {
625//
626 fMCGeo->Gsdvn(name, mother, ndiv, iaxis);
627}
628
629//______________________________________________________________________________
630void TFluka::Gsdvn2(const char *name, const char *mother, Int_t ndiv,
631 Int_t iaxis, Double_t c0i, Int_t numed) {
632//
633 fMCGeo->Gsdvn2(name, mother, ndiv, iaxis, c0i, numed);
634}
635
636//______________________________________________________________________________
637void TFluka::Gsdvt(const char *name, const char *mother, Double_t step,
638 Int_t iaxis, Int_t numed, Int_t ndvmx) {
639//
640 fMCGeo->Gsdvt(name, mother, step, iaxis, numed, ndvmx);
641}
642
643//______________________________________________________________________________
644void TFluka::Gsdvt2(const char *name, const char *mother, Double_t step,
645 Int_t iaxis, Double_t c0, Int_t numed, Int_t ndvmx) {
646//
647 fMCGeo->Gsdvt2(name, mother, step, iaxis, c0, numed, ndvmx);
648}
649
650//______________________________________________________________________________
651void TFluka::Gsord(const char * /*name*/, Int_t /*iax*/) {
652//
653// Nothing to do with TGeo
654}
655
656//______________________________________________________________________________
657void TFluka::Gspos(const char *name, Int_t nr, const char *mother,
658 Double_t x, Double_t y, Double_t z, Int_t irot,
659 const char *konly) {
660//
661 fMCGeo->Gspos(name, nr, mother, x, y, z, irot, konly);
662}
663
664//______________________________________________________________________________
665void TFluka::Gsposp(const char *name, Int_t nr, const char *mother,
666 Double_t x, Double_t y, Double_t z, Int_t irot,
667 const char *konly, Float_t *upar, Int_t np) {
668 //
669 fMCGeo->Gsposp(name, nr, mother, x, y, z, irot, konly, upar, np);
670}
671
672//______________________________________________________________________________
673void TFluka::Gsposp(const char *name, Int_t nr, const char *mother,
674 Double_t x, Double_t y, Double_t z, Int_t irot,
675 const char *konly, Double_t *upar, Int_t np) {
676 //
677 fMCGeo->Gsposp(name, nr, mother, x, y, z, irot, konly, upar, np);
678}
679
680//______________________________________________________________________________
681void TFluka::Gsbool(const char* /*onlyVolName*/, const char* /*manyVolName*/) {
682//
683// Nothing to do with TGeo
684}
685
686//______________________________________________________________________________
687void TFluka::SetCerenkov(Int_t itmed, Int_t npckov, Float_t* ppckov,
688 Float_t* absco, Float_t* effic, Float_t* rindex) {
689//
690// Set Cerenkov properties for medium itmed
691//
692// npckov: number of sampling points
693// ppckov: energy values
694// absco: absorption length
695// effic: quantum efficiency
696// rindex: refraction index
697//
698//
699//
700// Create object holding Cerenkov properties
701//
702 TFlukaCerenkov* cerenkovProperties = new TFlukaCerenkov(npckov, ppckov, absco, effic, rindex);
703//
704// Pass object to medium
705 TGeoMedium* medium = gGeoManager->GetMedium(itmed);
706 medium->SetCerenkovProperties(cerenkovProperties);
707}
708
709//______________________________________________________________________________
710void TFluka::SetCerenkov(Int_t /*itmed*/, Int_t /*npckov*/, Double_t * /*ppckov*/,
711 Double_t * /*absco*/, Double_t * /*effic*/, Double_t * /*rindex*/) {
712//
713// Not implemented with TGeo - what G4 did ? Any FLUKA card generated?
714 Warning("SetCerenkov", "Not implemented with TGeo");
715}
716
717// Euclid
718//______________________________________________________________________________
719void TFluka::WriteEuclid(const char* /*fileName*/, const char* /*topVol*/,
720 Int_t /*number*/, Int_t /*nlevel*/) {
721//
722// Not with TGeo
723 Warning("WriteEuclid", "Not implemented with TGeo");
724}
725
726
727
728//_____________________________________________________________________________
729// methods needed by the stepping
730//____________________________________________________________________________
731
732Int_t TFluka::GetMedium() const {
733//
734// Get the medium number for the current fluka region
735//
736 return fGeom->GetMedium(); // this I need to check due to remapping !!!
737}
738
739
740
741//____________________________________________________________________________
742// particle table usage
743// ID <--> PDG transformations
744//_____________________________________________________________________________
745Int_t TFluka::IdFromPDG(Int_t pdg) const
746{
747 //
748 // Return Fluka code from PDG and pseudo ENDF code
749
750 // Catch the feedback photons
751 if (pdg == 50000051) return (-1);
752 // MCIHAD() goes from pdg to fluka internal.
753 Int_t intfluka = mcihad(pdg);
754 // KPTOIP array goes from internal to official
755 return GetFlukaKPTOIP(intfluka);
756}
757
758//______________________________________________________________________________
759Int_t TFluka::PDGFromId(Int_t id) const
760{
761 //
762 // Return PDG code and pseudo ENDF code from Fluka code
f926898e 763 // Alpha He3 Triton Deuteron gen. ion opt. photon
764 Int_t idSpecial[6] = {10020040, 10020030, 10010030, 10010020, 10000000, 50000050};
829fb838 765 // IPTOKP array goes from official to internal
766
767 if (id == -1) {
768// Cerenkov photon
bd3d5c8a 769 if (fVerbosityLevel >= 3)
829fb838 770 printf("\n PDGFromId: Cerenkov Photon \n");
771 return 50000050;
772 }
773// Error id
774 if (id == 0 || id < -6 || id > 250) {
f926898e 775 if (fVerbosityLevel >= 3)
829fb838 776 printf("PDGFromId: Error id = 0\n");
777 return -1;
778 }
779// Good id
f926898e 780 if (id > 0) {
781 Int_t intfluka = GetFlukaIPTOKP(id);
782 if (intfluka == 0) {
783 if (fVerbosityLevel >= 3)
784 printf("PDGFromId: Error intfluka = 0: %d\n", id);
785 return -1;
786 } else if (intfluka < 0) {
787 if (fVerbosityLevel >= 3)
788 printf("PDGFromId: Error intfluka < 0: %d\n", id);
789 return -1;
790 }
791 if (fVerbosityLevel >= 3)
792 printf("mpdgha called with %d %d \n", id, intfluka);
793 // MPDGHA() goes from fluka internal to pdg.
794 return mpdgha(intfluka);
795 } else {
796 // ions and optical photons
797 return idSpecial[id + 6];
829fb838 798 }
829fb838 799}
800
bd3d5c8a 801void TFluka::StopTrack()
802{
803 // Set stopping conditions
804 // Works for photons and charged particles
805 fStopped = kTRUE;
806}
807
829fb838 808//_____________________________________________________________________________
809// methods for physics management
810//____________________________________________________________________________
811//
812// set methods
813//
814
1df5fa54 815void TFluka::SetProcess(const char* flagName, Int_t flagValue, Int_t imed)
829fb838 816{
817// Set process user flag for material imat
818//
1df5fa54 819//
820// Update if already in the list
821//
fb2cbbec 822 TIter next(fUserConfig);
1df5fa54 823 TFlukaConfigOption* proc;
824 while((proc = (TFlukaConfigOption*)next()))
825 {
fb2cbbec 826 if (proc->Medium() == imed) {
827 proc->SetProcess(flagName, flagValue);
828 return;
829 }
1df5fa54 830 }
fb2cbbec 831 proc = new TFlukaConfigOption(imed);
832 proc->SetProcess(flagName, flagValue);
833 fUserConfig->Add(proc);
834}
835
836//______________________________________________________________________________
837Bool_t TFluka::SetProcess(const char* flagName, Int_t flagValue)
838{
839// Set process user flag
1df5fa54 840//
1df5fa54 841//
fb2cbbec 842 SetProcess(flagName, flagValue, -1);
1df5fa54 843 return kTRUE;
829fb838 844}
845
846//______________________________________________________________________________
847void TFluka::SetCut(const char* cutName, Double_t cutValue, Int_t imed)
848{
849// Set user cut value for material imed
850//
fb2cbbec 851 TIter next(fUserConfig);
852 TFlukaConfigOption* proc;
853 while((proc = (TFlukaConfigOption*)next()))
854 {
855 if (proc->Medium() == imed) {
856 proc->SetCut(cutName, cutValue);
857 return;
858 }
859 }
860
861 proc = new TFlukaConfigOption(imed);
862 proc->SetCut(cutName, cutValue);
863 fUserConfig->Add(proc);
829fb838 864}
865
866//______________________________________________________________________________
867Bool_t TFluka::SetCut(const char* cutName, Double_t cutValue)
868{
869// Set user cut value
870//
1df5fa54 871//
fb2cbbec 872 SetCut(cutName, cutValue, -1);
873 return kTRUE;
829fb838 874}
875
b496f27c 876void TFluka::SetUserScoring(const char* option, Int_t npar, Float_t what[12])
877{
878//
879// Ads a user scoring option to th list
880//
881 TFlukaScoringOption* opt = new TFlukaScoringOption(option, "User Scoring", npar, what);
882 fUserScore->Add(opt);
883}
884
885
829fb838 886//______________________________________________________________________________
887Double_t TFluka::Xsec(char*, Double_t, Int_t, Int_t)
888{
889 printf("WARNING: Xsec not yet implemented !\n"); return -1.;
890}
891
892
893//______________________________________________________________________________
894void TFluka::InitPhysics()
895{
896//
897// Physics initialisation with preparation of FLUKA input cards
898//
fb2cbbec 899 printf("=>InitPhysics\n");
829fb838 900
fb2cbbec 901// Construct file names
902 FILE *pFlukaVmcCoreInp, *pFlukaVmcFlukaMat, *pFlukaVmcInp;
903 TString sFlukaVmcCoreInp = getenv("ALICE_ROOT");
904 sFlukaVmcCoreInp +="/TFluka/input/";
905 TString sFlukaVmcTmp = "flukaMat.inp";
906 TString sFlukaVmcInp = GetInputFileName();
907 sFlukaVmcCoreInp += GetCoreInputFileName();
908
909// Open files
910 if ((pFlukaVmcCoreInp = fopen(sFlukaVmcCoreInp.Data(),"r")) == NULL) {
911 printf("\nCannot open file %s\n",sFlukaVmcCoreInp.Data());
912 exit(1);
913 }
914 if ((pFlukaVmcFlukaMat = fopen(sFlukaVmcTmp.Data(),"r")) == NULL) {
915 printf("\nCannot open file %s\n",sFlukaVmcTmp.Data());
916 exit(1);
917 }
918 if ((pFlukaVmcInp = fopen(sFlukaVmcInp.Data(),"w")) == NULL) {
919 printf("\nCannot open file %s\n",sFlukaVmcInp.Data());
920 exit(1);
921 }
829fb838 922
fb2cbbec 923// Copy core input file
924 Char_t sLine[255];
925 Float_t fEventsPerRun;
829fb838 926
fb2cbbec 927 while ((fgets(sLine,255,pFlukaVmcCoreInp)) != NULL) {
928 if (strncmp(sLine,"GEOEND",6) != 0)
929 fprintf(pFlukaVmcInp,"%s",sLine); // copy until GEOEND card
930 else {
931 fprintf(pFlukaVmcInp,"GEOEND\n"); // add GEOEND card
932 goto flukamat;
829fb838 933 }
fb2cbbec 934 } // end of while until GEOEND card
935
829fb838 936
fb2cbbec 937 flukamat:
938 while ((fgets(sLine,255,pFlukaVmcFlukaMat)) != NULL) { // copy flukaMat.inp file
939 fprintf(pFlukaVmcInp,"%s\n",sLine);
940 }
941
942 while ((fgets(sLine,255,pFlukaVmcCoreInp)) != NULL) {
943 if (strncmp(sLine,"START",5) != 0)
944 fprintf(pFlukaVmcInp,"%s\n",sLine);
945 else {
946 sscanf(sLine+10,"%10f",&fEventsPerRun);
947 goto fin;
948 }
949 } //end of while until START card
950
951 fin:
829fb838 952
fb2cbbec 953 // Pass information to configuration objects
829fb838 954
fb2cbbec 955 Float_t fLastMaterial = fGeom->GetLastMaterialIndex();
956 TFlukaConfigOption::SetStaticInfo(pFlukaVmcInp, 3, fLastMaterial, fGeom);
957
958 TIter next(fUserConfig);
959 TFlukaConfigOption* proc;
960 while((proc = (TFlukaConfigOption*)next())) proc->WriteFlukaInputCards();
961
829fb838 962// Add START and STOP card
3b8c325d 963 fprintf(pFlukaVmcInp,"START %10.1f\n",fEventsPerRun);
964 fprintf(pFlukaVmcInp,"STOP \n");
829fb838 965
966
967// Close files
3b8c325d 968 fclose(pFlukaVmcCoreInp);
969 fclose(pFlukaVmcFlukaMat);
970 fclose(pFlukaVmcInp);
fb2cbbec 971
972
973//
974// Initialisation needed for Cerenkov photon production and transport
975 TObjArray *matList = GetFlukaMaterials();
976 Int_t nmaterial = matList->GetEntriesFast();
977 fMaterials = new Int_t[nmaterial+3];
978
979 for (Int_t im = 0; im < nmaterial; im++)
980 {
981 TGeoMaterial* material = dynamic_cast<TGeoMaterial*> (matList->At(im));
982 Int_t idmat = material->GetIndex();
983 fMaterials[idmat] = im;
984 }
829fb838 985} // end of InitPhysics
986
987
988//______________________________________________________________________________
07f5b33e 989void TFluka::SetMaxStep(Double_t step)
829fb838 990{
07f5b33e 991// Set the maximum step size
992 if (step > 1.e4) return;
993
994 Int_t mreg, latt;
995 fGeom->GetCurrentRegion(mreg, latt);
996 STEPSZ.stepmx[mreg - 1] = step;
829fb838 997}
998
2f09b80e 999
1000Double_t TFluka::MaxStep() const
1001{
1002// Return the maximum for current medium
1003 Int_t mreg, latt;
1004 fGeom->GetCurrentRegion(mreg, latt);
1005 return (STEPSZ.stepmx[mreg - 1]);
1006}
1007
829fb838 1008//______________________________________________________________________________
1009void TFluka::SetMaxNStep(Int_t)
1010{
1011// SetMaxNStep is dummy procedure in TFluka !
1012 if (fVerbosityLevel >=3)
1013 cout << "SetMaxNStep is dummy procedure in TFluka !" << endl;
1014}
1015
1016//______________________________________________________________________________
1017void TFluka::SetUserDecay(Int_t)
1018{
1019// SetUserDecay is dummy procedure in TFluka !
1020 if (fVerbosityLevel >=3)
1021 cout << "SetUserDecay is dummy procedure in TFluka !" << endl;
1022}
1023
1024//
1025// dynamic properties
1026//
1027//______________________________________________________________________________
1028void TFluka::TrackPosition(TLorentzVector& position) const
1029{
1030// Return the current position in the master reference frame of the
1031// track being transported
1032// TRACKR.atrack = age of the particle
1033// TRACKR.xtrack = x-position of the last point
1034// TRACKR.ytrack = y-position of the last point
1035// TRACKR.ztrack = z-position of the last point
1036 Int_t caller = GetCaller();
1037 if (caller == 3 || caller == 6 || caller == 11 || caller == 12) { //bxdraw,endraw,usdraw
1038 position.SetX(GetXsco());
1039 position.SetY(GetYsco());
1040 position.SetZ(GetZsco());
1041 position.SetT(TRACKR.atrack);
1042 }
1043 else if (caller == 4) { // mgdraw
1044 position.SetX(TRACKR.xtrack[TRACKR.ntrack]);
1045 position.SetY(TRACKR.ytrack[TRACKR.ntrack]);
1046 position.SetZ(TRACKR.ztrack[TRACKR.ntrack]);
1047 position.SetT(TRACKR.atrack);
1048 }
1049 else if (caller == 5) { // sodraw
1050 position.SetX(TRACKR.xtrack[TRACKR.ntrack]);
1051 position.SetY(TRACKR.ytrack[TRACKR.ntrack]);
1052 position.SetZ(TRACKR.ztrack[TRACKR.ntrack]);
1053 position.SetT(0);
1054 }
1055 else
1056 Warning("TrackPosition","position not available");
1057}
1058
1059//______________________________________________________________________________
1060void TFluka::TrackPosition(Double_t& x, Double_t& y, Double_t& z) const
1061{
1062// Return the current position in the master reference frame of the
1063// track being transported
1064// TRACKR.atrack = age of the particle
1065// TRACKR.xtrack = x-position of the last point
1066// TRACKR.ytrack = y-position of the last point
1067// TRACKR.ztrack = z-position of the last point
1068 Int_t caller = GetCaller();
1069 if (caller == 3 || caller == 6 || caller == 11 || caller == 12) { //bxdraw,endraw,usdraw
1070 x = GetXsco();
1071 y = GetYsco();
1072 z = GetZsco();
1073 }
1074 else if (caller == 4 || caller == 5) { // mgdraw, sodraw
1075 x = TRACKR.xtrack[TRACKR.ntrack];
1076 y = TRACKR.ytrack[TRACKR.ntrack];
1077 z = TRACKR.ztrack[TRACKR.ntrack];
1078 }
1079 else
1080 Warning("TrackPosition","position not available");
1081}
1082
1083//______________________________________________________________________________
1084void TFluka::TrackMomentum(TLorentzVector& momentum) const
1085{
1086// Return the direction and the momentum (GeV/c) of the track
1087// currently being transported
1088// TRACKR.ptrack = momentum of the particle (not always defined, if
1089// < 0 must be obtained from etrack)
1090// TRACKR.cx,y,ztrck = direction cosines of the current particle
1091// TRACKR.etrack = total energy of the particle
1092// TRACKR.jtrack = identity number of the particle
1093// PAPROP.am[TRACKR.jtrack] = particle mass in gev
1094 Int_t caller = GetCaller();
1095 if (caller != 2) { // not eedraw
1096 if (TRACKR.ptrack >= 0) {
1097 momentum.SetPx(TRACKR.ptrack*TRACKR.cxtrck);
1098 momentum.SetPy(TRACKR.ptrack*TRACKR.cytrck);
1099 momentum.SetPz(TRACKR.ptrack*TRACKR.cztrck);
1100 momentum.SetE(TRACKR.etrack);
1101 return;
1102 }
1103 else {
1104 Double_t p = sqrt(TRACKR.etrack*TRACKR.etrack - PAPROP.am[TRACKR.jtrack+6]*PAPROP.am[TRACKR.jtrack+6]);
1105 momentum.SetPx(p*TRACKR.cxtrck);
1106 momentum.SetPy(p*TRACKR.cytrck);
1107 momentum.SetPz(p*TRACKR.cztrck);
1108 momentum.SetE(TRACKR.etrack);
1109 return;
1110 }
1111 }
1112 else
1113 Warning("TrackMomentum","momentum not available");
1114}
1115
1116//______________________________________________________________________________
1117void TFluka::TrackMomentum(Double_t& px, Double_t& py, Double_t& pz, Double_t& e) const
1118{
1119// Return the direction and the momentum (GeV/c) of the track
1120// currently being transported
1121// TRACKR.ptrack = momentum of the particle (not always defined, if
1122// < 0 must be obtained from etrack)
1123// TRACKR.cx,y,ztrck = direction cosines of the current particle
1124// TRACKR.etrack = total energy of the particle
1125// TRACKR.jtrack = identity number of the particle
1126// PAPROP.am[TRACKR.jtrack] = particle mass in gev
1127 Int_t caller = GetCaller();
1128 if (caller != 2) { // not eedraw
1129 if (TRACKR.ptrack >= 0) {
1130 px = TRACKR.ptrack*TRACKR.cxtrck;
1131 py = TRACKR.ptrack*TRACKR.cytrck;
1132 pz = TRACKR.ptrack*TRACKR.cztrck;
1133 e = TRACKR.etrack;
1134 return;
1135 }
1136 else {
1137 Double_t p = sqrt(TRACKR.etrack*TRACKR.etrack - PAPROP.am[TRACKR.jtrack+6]*PAPROP.am[TRACKR.jtrack+6]);
1138 px = p*TRACKR.cxtrck;
1139 py = p*TRACKR.cytrck;
1140 pz = p*TRACKR.cztrck;
1141 e = TRACKR.etrack;
1142 return;
1143 }
1144 }
1145 else
1146 Warning("TrackMomentum","momentum not available");
1147}
1148
1149//______________________________________________________________________________
1150Double_t TFluka::TrackStep() const
1151{
1152// Return the length in centimeters of the current step
1153// TRACKR.ctrack = total curved path
1154 Int_t caller = GetCaller();
1155 if (caller == 11 || caller==12 || caller == 3 || caller == 6) //bxdraw,endraw,usdraw
1156 return 0.0;
1157 else if (caller == 4) //mgdraw
1158 return TRACKR.ctrack;
1159 else
1160 return -1.0;
1161}
1162
1163//______________________________________________________________________________
1164Double_t TFluka::TrackLength() const
1165{
1166// TRACKR.cmtrck = cumulative curved path since particle birth
1167 Int_t caller = GetCaller();
1168 if (caller == 11 || caller==12 || caller == 3 || caller == 4 || caller == 6) //bxdraw,endraw,mgdraw,usdraw
1169 return TRACKR.cmtrck;
1170 else
1171 return -1.0;
1172}
1173
1174//______________________________________________________________________________
1175Double_t TFluka::TrackTime() const
1176{
1177// Return the current time of flight of the track being transported
1178// TRACKR.atrack = age of the particle
1179 Int_t caller = GetCaller();
1180 if (caller == 11 || caller==12 || caller == 3 || caller == 4 || caller == 6) //bxdraw,endraw,mgdraw,usdraw
1181 return TRACKR.atrack;
1182 else
1183 return -1;
1184}
1185
1186//______________________________________________________________________________
1187Double_t TFluka::Edep() const
1188{
1189// Energy deposition
1190// if TRACKR.ntrack = 0, TRACKR.mtrack = 0:
1191// -->local energy deposition (the value and the point are not recorded in TRACKR)
1192// but in the variable "rull" of the procedure "endraw.cxx"
1193// if TRACKR.ntrack > 0, TRACKR.mtrack = 0:
1194// -->no energy loss along the track
1195// if TRACKR.ntrack > 0, TRACKR.mtrack > 0:
1196// -->energy loss distributed along the track
07f5b33e 1197// TRACKR.dtrack = energy deposition of the jth deposition event
829fb838 1198
1199 // If coming from bxdraw we have 2 steps of 0 length and 0 edep
1200 Int_t caller = GetCaller();
1201 if (caller == 11 || caller==12) return 0.0;
1202 Double_t sum = 0;
1203 for ( Int_t j=0;j<TRACKR.mtrack;j++) {
1204 sum +=TRACKR.dtrack[j];
1205 }
1206 if (TRACKR.ntrack == 0 && TRACKR.mtrack == 0)
1207 return fRull + sum;
1208 else {
1209 return sum;
1210 }
1211}
1212
1213//______________________________________________________________________________
1214Int_t TFluka::TrackPid() const
1215{
1216// Return the id of the particle transported
1217// TRACKR.jtrack = identity number of the particle
1218 Int_t caller = GetCaller();
f926898e 1219 if (caller != 2) { // not eedraw
1220 return PDGFromId(TRACKR.jtrack);
1221 }
829fb838 1222 else
1223 return -1000;
1224}
1225
1226//______________________________________________________________________________
1227Double_t TFluka::TrackCharge() const
1228{
1229// Return charge of the track currently transported
1230// PAPROP.ichrge = electric charge of the particle
1231// TRACKR.jtrack = identity number of the particle
1232 Int_t caller = GetCaller();
1233 if (caller != 2) // not eedraw
1234 return PAPROP.ichrge[TRACKR.jtrack+6];
1235 else
1236 return -1000.0;
1237}
1238
1239//______________________________________________________________________________
1240Double_t TFluka::TrackMass() const
1241{
1242// PAPROP.am = particle mass in GeV
1243// TRACKR.jtrack = identity number of the particle
1244 Int_t caller = GetCaller();
1245 if (caller != 2) // not eedraw
1246 return PAPROP.am[TRACKR.jtrack+6];
1247 else
1248 return -1000.0;
1249}
1250
1251//______________________________________________________________________________
1252Double_t TFluka::Etot() const
1253{
1254// TRACKR.etrack = total energy of the particle
1255 Int_t caller = GetCaller();
1256 if (caller != 2) // not eedraw
1257 return TRACKR.etrack;
1258 else
1259 return -1000.0;
1260}
1261
1262//
1263// track status
1264//
1265//______________________________________________________________________________
1266Bool_t TFluka::IsNewTrack() const
1267{
1268// Return true for the first call of Stepping()
1269 return fTrackIsNew;
1270}
1271
0dabe425 1272void TFluka::SetTrackIsNew(Bool_t flag)
1273{
1274// Return true for the first call of Stepping()
1275 fTrackIsNew = flag;
1276
1277}
1278
1279
829fb838 1280//______________________________________________________________________________
1281Bool_t TFluka::IsTrackInside() const
1282{
1283// True if the track is not at the boundary of the current volume
1284// In Fluka a step is always inside one kind of material
1285// If the step would go behind the region of one material,
1286// it will be shortened to reach only the boundary.
1287// Therefore IsTrackInside() is always true.
1288 Int_t caller = GetCaller();
1289 if (caller == 11 || caller==12) // bxdraw
1290 return 0;
1291 else
1292 return 1;
1293}
1294
1295//______________________________________________________________________________
1296Bool_t TFluka::IsTrackEntering() const
1297{
1298// True if this is the first step of the track in the current volume
1299
1300 Int_t caller = GetCaller();
1301 if (caller == 11) // bxdraw entering
1302 return 1;
1303 else return 0;
1304}
1305
1306//______________________________________________________________________________
1307Bool_t TFluka::IsTrackExiting() const
1308{
1309// True if track is exiting volume
1310//
1311 Int_t caller = GetCaller();
1312 if (caller == 12) // bxdraw exiting
1313 return 1;
1314 else return 0;
1315}
1316
1317//______________________________________________________________________________
1318Bool_t TFluka::IsTrackOut() const
1319{
1320// True if the track is out of the setup
1321// means escape
1322// Icode = 14: escape - call from Kaskad
1323// Icode = 23: escape - call from Emfsco
1324// Icode = 32: escape - call from Kasneu
1325// Icode = 40: escape - call from Kashea
1326// Icode = 51: escape - call from Kasoph
1327 if (fIcode == 14 ||
1328 fIcode == 23 ||
1329 fIcode == 32 ||
1330 fIcode == 40 ||
1331 fIcode == 51) return 1;
1332 else return 0;
1333}
1334
1335//______________________________________________________________________________
1336Bool_t TFluka::IsTrackDisappeared() const
1337{
1338// means all inelastic interactions and decays
1339// fIcode from usdraw
1340 if (fIcode == 101 || // inelastic interaction
1341 fIcode == 102 || // particle decay
0dabe425 1342 fIcode == 103 || // delta ray generation by hadron
1343 fIcode == 104 || // direct pair production
1344 fIcode == 105 || // bremsstrahlung (muon)
1345 fIcode == 208 || // bremsstrahlung (electron)
829fb838 1346 fIcode == 214 || // in-flight annihilation
1347 fIcode == 215 || // annihilation at rest
1348 fIcode == 217 || // pair production
0dabe425 1349 fIcode == 219 || // Compton scattering
1350 fIcode == 221 || // Photoelectric effect
1351 fIcode == 300 || // hadronic interaction
1352 fIcode == 400 // delta-ray
1353 ) return 1;
829fb838 1354 else return 0;
1355}
1356
1357//______________________________________________________________________________
1358Bool_t TFluka::IsTrackStop() const
1359{
1360// True if the track energy has fallen below the threshold
1361// means stopped by signal or below energy threshold
1362// Icode = 12: stopping particle - call from Kaskad
1363// Icode = 15: time kill - call from Kaskad
1364// Icode = 21: below threshold, iarg=1 - call from Emfsco
1365// Icode = 22: below threshold, iarg=2 - call from Emfsco
1366// Icode = 24: time kill - call from Emfsco
1367// Icode = 31: below threshold - call from Kasneu
1368// Icode = 33: time kill - call from Kasneu
1369// Icode = 41: time kill - call from Kashea
1370// Icode = 52: time kill - call from Kasoph
1371 if (fIcode == 12 ||
1372 fIcode == 15 ||
1373 fIcode == 21 ||
1374 fIcode == 22 ||
1375 fIcode == 24 ||
1376 fIcode == 31 ||
1377 fIcode == 33 ||
1378 fIcode == 41 ||
1379 fIcode == 52) return 1;
1380 else return 0;
1381}
1382
1383//______________________________________________________________________________
1384Bool_t TFluka::IsTrackAlive() const
1385{
1386// means not disappeared or not out
1387 if (IsTrackDisappeared() || IsTrackOut() ) return 0;
1388 else return 1;
1389}
1390
1391//
1392// secondaries
1393//
1394
1395//______________________________________________________________________________
1396Int_t TFluka::NSecondaries() const
1397
1398{
1399// Number of secondary particles generated in the current step
1400// FINUC.np = number of secondaries except light and heavy ions
1401// FHEAVY.npheav = number of secondaries for light and heavy secondary ions
7b203b6e 1402 Int_t caller = GetCaller();
1403 if (caller == 6) // valid only after usdraw
1404 return FINUC.np + FHEAVY.npheav;
1405 else if (caller == 50) {
1406 // Cerenkov Photon production
1407 return fNCerenkov;
1408 }
829fb838 1409 return 0;
1410} // end of NSecondaries
1411
1412//______________________________________________________________________________
1413void TFluka::GetSecondary(Int_t isec, Int_t& particleId,
1414 TLorentzVector& position, TLorentzVector& momentum)
1415{
1416// Copy particles from secondary stack to vmc stack
1417//
1418
7b203b6e 1419 Int_t caller = GetCaller();
1420 if (caller == 6) { // valid only after usdraw
1421 if (FINUC.np > 0) {
1422 // Hadronic interaction
1423 if (isec >= 0 && isec < FINUC.np) {
1424 particleId = PDGFromId(FINUC.kpart[isec]);
1425 position.SetX(fXsco);
1426 position.SetY(fYsco);
1427 position.SetZ(fZsco);
1428 position.SetT(TRACKR.atrack);
1429 momentum.SetPx(FINUC.plr[isec]*FINUC.cxr[isec]);
1430 momentum.SetPy(FINUC.plr[isec]*FINUC.cyr[isec]);
1431 momentum.SetPz(FINUC.plr[isec]*FINUC.czr[isec]);
1432 momentum.SetE(FINUC.tki[isec] + PAPROP.am[FINUC.kpart[isec]+6]);
1433 }
1434 else if (isec >= FINUC.np && isec < FINUC.np + FHEAVY.npheav) {
1435 Int_t jsec = isec - FINUC.np;
1436 particleId = FHEAVY.kheavy[jsec]; // this is Fluka id !!!
1437 position.SetX(fXsco);
1438 position.SetY(fYsco);
1439 position.SetZ(fZsco);
1440 position.SetT(TRACKR.atrack);
1441 momentum.SetPx(FHEAVY.pheavy[jsec]*FHEAVY.cxheav[jsec]);
1442 momentum.SetPy(FHEAVY.pheavy[jsec]*FHEAVY.cyheav[jsec]);
1443 momentum.SetPz(FHEAVY.pheavy[jsec]*FHEAVY.czheav[jsec]);
1444 if (FHEAVY.tkheav[jsec] >= 3 && FHEAVY.tkheav[jsec] <= 6)
1445 momentum.SetE(FHEAVY.tkheav[jsec] + PAPROP.am[jsec+6]);
1446 else if (FHEAVY.tkheav[jsec] > 6)
1447 momentum.SetE(FHEAVY.tkheav[jsec] + FHEAVY.amnhea[jsec]); // to be checked !!!
1448 }
1449 else
1450 Warning("GetSecondary","isec out of range");
1451 }
1452 } else if (caller == 50) {
1453 Int_t index = OPPHST.lstopp - isec;
1454 position.SetX(OPPHST.xoptph[index]);
1455 position.SetY(OPPHST.yoptph[index]);
1456 position.SetZ(OPPHST.zoptph[index]);
1457 position.SetT(OPPHST.agopph[index]);
1458 Double_t p = OPPHST.poptph[index];
1459
1460 momentum.SetPx(p * OPPHST.txopph[index]);
1461 momentum.SetPy(p * OPPHST.tyopph[index]);
1462 momentum.SetPz(p * OPPHST.tzopph[index]);
1463 momentum.SetE(p);
829fb838 1464 }
1465 else
7b203b6e 1466 Warning("GetSecondary","no secondaries available");
1467
829fb838 1468} // end of GetSecondary
1469
7b203b6e 1470
829fb838 1471//______________________________________________________________________________
1472TMCProcess TFluka::ProdProcess(Int_t) const
1473
1474{
1475// Name of the process that has produced the secondary particles
1476// in the current step
0dabe425 1477
1478 Int_t mugamma = (TRACKR.jtrack == 7 || TRACKR.jtrack == 10 || TRACKR.jtrack == 11);
1479
b496f27c 1480 if (fIcode == 102) return kPDecay;
0dabe425 1481 else if (fIcode == 104 || fIcode == 217) return kPPair;
b496f27c 1482 else if (fIcode == 219) return kPCompton;
1483 else if (fIcode == 221) return kPPhotoelectric;
0dabe425 1484 else if (fIcode == 105 || fIcode == 208) return kPBrem;
1485 else if (fIcode == 103 || fIcode == 400) return kPDeltaRay;
1486 else if (fIcode == 210 || fIcode == 212) return kPDeltaRay;
1487 else if (fIcode == 214 || fIcode == 215) return kPAnnihilation;
b496f27c 1488 else if (fIcode == 101) return kPHadronic;
829fb838 1489 else if (fIcode == 101) {
b496f27c 1490 if (!mugamma) return kPHadronic;
1491 else if (TRACKR.jtrack == 7) return kPPhotoFission;
1492 else return kPMuonNuclear;
829fb838 1493 }
b496f27c 1494 else if (fIcode == 225) return kPRayleigh;
829fb838 1495// Fluka codes 100, 300 and 400 still to be investigasted
b496f27c 1496 else return kPNoProcess;
829fb838 1497}
1498
829fb838 1499
b496f27c 1500Int_t TFluka::StepProcesses(TArrayI &proc) const
1501{
1502 //
1503 // Return processes active in the current step
1504 //
1505 proc.Set(1);
1506 TMCProcess iproc;
1507 switch (fIcode) {
1508 case 15:
1509 case 24:
1510 case 33:
1511 case 41:
1512 case 52:
1513 iproc = kPTOFlimit;
1514 break;
1515 case 12:
1516 case 14:
1517 case 21:
1518 case 22:
1519 case 23:
1520 case 31:
1521 case 32:
1522 case 40:
1523 case 51:
6fd5baa4 1524 iproc = kPStop;
b496f27c 1525 break;
1526 case 50:
1527 iproc = kPLightAbsorption;
1528 break;
6fd5baa4 1529 case 59:
1530 iproc = kPLightRefraction;
b496f27c 1531 case 20:
1532 iproc = kPPhotoelectric;
1533 break;
1534 default:
1535 iproc = ProdProcess(0);
1536 }
07f5b33e 1537 proc[0] = iproc;
b496f27c 1538 return 1;
1539}
829fb838 1540//______________________________________________________________________________
1541Int_t TFluka::VolId2Mate(Int_t id) const
1542{
1543//
1544// Returns the material number for a given volume ID
1545//
1546 return fMCGeo->VolId2Mate(id);
1547}
1548
1549//______________________________________________________________________________
1550const char* TFluka::VolName(Int_t id) const
1551{
1552//
1553// Returns the volume name for a given volume ID
1554//
1555 return fMCGeo->VolName(id);
1556}
1557
1558//______________________________________________________________________________
1559Int_t TFluka::VolId(const Text_t* volName) const
1560{
1561//
1562// Converts from volume name to volume ID.
1563// Time consuming. (Only used during set-up)
1564// Could be replaced by hash-table
1565//
09cd6497 1566 char sname[20];
1567 Int_t len;
1568 strncpy(sname, volName, len = strlen(volName));
1569 sname[len] = 0;
1570 while (sname[len - 1] == ' ') sname[--len] = 0;
1571 return fMCGeo->VolId(sname);
829fb838 1572}
1573
1574//______________________________________________________________________________
1575Int_t TFluka::CurrentVolID(Int_t& copyNo) const
1576{
1577//
1578// Return the logical id and copy number corresponding to the current fluka region
1579//
1580 if (gGeoManager->IsOutside()) return 0;
1581 TGeoNode *node = gGeoManager->GetCurrentNode();
1582 copyNo = node->GetNumber();
1583 Int_t id = node->GetVolume()->GetNumber();
1584 return id;
1585}
1586
1587//______________________________________________________________________________
1588Int_t TFluka::CurrentVolOffID(Int_t off, Int_t& copyNo) const
1589{
1590//
1591// Return the logical id and copy number of off'th mother
1592// corresponding to the current fluka region
1593//
1594 if (off<0 || off>gGeoManager->GetLevel()) return 0;
1595 if (off==0) return CurrentVolID(copyNo);
1596 TGeoNode *node = gGeoManager->GetMother(off);
1597 if (!node) return 0;
1598 copyNo = node->GetNumber();
1599 return node->GetVolume()->GetNumber();
1600}
1601
1602//______________________________________________________________________________
1603const char* TFluka::CurrentVolName() const
1604{
1605//
1606// Return the current volume name
1607//
1608 if (gGeoManager->IsOutside()) return 0;
1609 return gGeoManager->GetCurrentVolume()->GetName();
1610}
1611
1612//______________________________________________________________________________
1613const char* TFluka::CurrentVolOffName(Int_t off) const
1614{
1615//
1616// Return the volume name of the off'th mother of the current volume
1617//
1618 if (off<0 || off>gGeoManager->GetLevel()) return 0;
1619 if (off==0) return CurrentVolName();
1620 TGeoNode *node = gGeoManager->GetMother(off);
1621 if (!node) return 0;
1622 return node->GetVolume()->GetName();
1623}
1624
1625//______________________________________________________________________________
a60813de 1626Int_t TFluka::CurrentMaterial(Float_t & a, Float_t & z,
1627 Float_t & dens, Float_t & radl, Float_t & absl) const
829fb838 1628{
1629//
a60813de 1630// Return the current medium number and material properties
829fb838 1631//
1632 Int_t copy;
1633 Int_t id = TFluka::CurrentVolID(copy);
1634 Int_t med = TFluka::VolId2Mate(id);
a60813de 1635 TGeoVolume* vol = gGeoManager->GetCurrentVolume();
1636 TGeoMaterial* mat = vol->GetMaterial();
1637 a = mat->GetA();
1638 z = mat->GetZ();
1639 dens = mat->GetDensity();
1640 radl = mat->GetRadLen();
1641 absl = mat->GetIntLen();
1642
829fb838 1643 return med;
1644}
1645
1646//______________________________________________________________________________
1647void TFluka::Gmtod(Float_t* xm, Float_t* xd, Int_t iflag)
1648{
1649// Transforms a position from the world reference frame
1650// to the current volume reference frame.
1651//
1652// Geant3 desription:
1653// ==================
1654// Computes coordinates XD (in DRS)
1655// from known coordinates XM in MRS
1656// The local reference system can be initialized by
1657// - the tracking routines and GMTOD used in GUSTEP
1658// - a call to GMEDIA(XM,NUMED)
1659// - a call to GLVOLU(NLEVEL,NAMES,NUMBER,IER)
1660// (inverse routine is GDTOM)
1661//
1662// If IFLAG=1 convert coordinates
1663// IFLAG=2 convert direction cosinus
1664//
1665// ---
1666 Double_t xmL[3], xdL[3];
1667 Int_t i;
1668 for (i=0;i<3;i++) xmL[i]=xm[i];
1669 if (iflag == 1) gGeoManager->MasterToLocal(xmL,xdL);
1670 else gGeoManager->MasterToLocalVect(xmL,xdL);
1671 for (i=0;i<3;i++) xd[i] = xdL[i];
1672}
1673
1674//______________________________________________________________________________
1675void TFluka::Gmtod(Double_t* xm, Double_t* xd, Int_t iflag)
1676{
1677 if (iflag == 1) gGeoManager->MasterToLocal(xm,xd);
1678 else gGeoManager->MasterToLocalVect(xm,xd);
1679}
1680
1681//______________________________________________________________________________
1682void TFluka::Gdtom(Float_t* xd, Float_t* xm, Int_t iflag)
1683{
1684// Transforms a position from the current volume reference frame
1685// to the world reference frame.
1686//
1687// Geant3 desription:
1688// ==================
1689// Computes coordinates XM (Master Reference System
1690// knowing the coordinates XD (Detector Ref System)
1691// The local reference system can be initialized by
1692// - the tracking routines and GDTOM used in GUSTEP
1693// - a call to GSCMED(NLEVEL,NAMES,NUMBER)
1694// (inverse routine is GMTOD)
1695//
1696// If IFLAG=1 convert coordinates
1697// IFLAG=2 convert direction cosinus
1698//
1699// ---
1700 Double_t xmL[3], xdL[3];
1701 Int_t i;
1702 for (i=0;i<3;i++) xdL[i] = xd[i];
1703 if (iflag == 1) gGeoManager->LocalToMaster(xdL,xmL);
1704 else gGeoManager->LocalToMasterVect(xdL,xmL);
1705 for (i=0;i<3;i++) xm[i]=xmL[i];
1706}
1707
1708//______________________________________________________________________________
1709void TFluka::Gdtom(Double_t* xd, Double_t* xm, Int_t iflag)
1710{
1711 if (iflag == 1) gGeoManager->LocalToMaster(xd,xm);
1712 else gGeoManager->LocalToMasterVect(xd,xm);
1713}
1714
1715//______________________________________________________________________________
1716TObjArray *TFluka::GetFlukaMaterials()
1717{
1718 return fGeom->GetMatList();
1719}
1720
1721//______________________________________________________________________________
1722void TFluka::SetMreg(Int_t l)
1723{
1724// Set current fluka region
1725 fCurrentFlukaRegion = l;
1726 fGeom->SetMreg(l);
1727}
1728
1729
b496f27c 1730
1731
1732TString TFluka::ParticleName(Int_t pdg) const
1733{
1734 // Return particle name for particle with pdg code pdg.
1735 Int_t ifluka = IdFromPDG(pdg);
1736 return TString((CHPPRP.btype[ifluka+6]), 8);
1737}
1738
1739
1740Double_t TFluka::ParticleMass(Int_t pdg) const
1741{
1742 // Return particle mass for particle with pdg code pdg.
1743 Int_t ifluka = IdFromPDG(pdg);
1744 return (PAPROP.am[ifluka+6]);
1745}
1746
1747Double_t TFluka::ParticleCharge(Int_t pdg) const
1748{
1749 // Return particle charge for particle with pdg code pdg.
1750 Int_t ifluka = IdFromPDG(pdg);
1751 return Double_t(PAPROP.ichrge[ifluka+6]);
1752}
1753
1754Double_t TFluka::ParticleLifeTime(Int_t pdg) const
1755{
1756 // Return particle lifetime for particle with pdg code pdg.
1757 Int_t ifluka = IdFromPDG(pdg);
1758 return (PAPROP.thalf[ifluka+6]);
1759}
1760
1761void TFluka::Gfpart(Int_t pdg, char* name, Int_t& type, Float_t& mass, Float_t& charge, Float_t& tlife)
1762{
1763 // Retrieve particle properties for particle with pdg code pdg.
1764
1765 strcpy(name, ParticleName(pdg).Data());
1766 type = ParticleMCType(pdg);
1767 mass = ParticleMass(pdg);
1768 charge = ParticleCharge(pdg);
1769 tlife = ParticleLifeTime(pdg);
1770}
1771
1772
1773
3a625972 1774#define pushcerenkovphoton pushcerenkovphoton_
7b203b6e 1775#define usersteppingckv usersteppingckv_
3a625972 1776
1777
1778extern "C" {
1779 void pushcerenkovphoton(Double_t & px, Double_t & py, Double_t & pz, Double_t & e,
1780 Double_t & vx, Double_t & vy, Double_t & vz, Double_t & tof,
1781 Double_t & polx, Double_t & poly, Double_t & polz, Double_t & wgt, Int_t& ntr)
1782 {
1783 //
1784 // Pushes one cerenkov photon to the stack
1785 //
1786
1787 TFluka* fluka = (TFluka*) gMC;
1788 TVirtualMCStack* cppstack = fluka->GetStack();
bd3d5c8a 1789 Int_t parent = TRACKR.ispusr[mkbmx2-1];
921e0994 1790 cppstack->PushTrack(0, parent, 50000050,
3a625972 1791 px, py, pz, e,
1792 vx, vy, vz, tof,
1793 polx, poly, polz,
1794 kPCerenkov, ntr, wgt, 0);
1795 }
7b203b6e 1796
1797 void usersteppingckv(Int_t & nphot, Int_t & mreg, Double_t & x, Double_t & y, Double_t & z)
1798 {
1799 //
1800 // Calls stepping in order to signal cerenkov production
1801 //
1802 TFluka *fluka = (TFluka*)gMC;
1803 fluka->SetMreg(mreg);
1804 fluka->SetXsco(x);
1805 fluka->SetYsco(y);
1806 fluka->SetZsco(z);
1807 fluka->SetNCerenkov(nphot);
1808 fluka->SetCaller(50);
1809 printf("userstepping ckv: %10d %10d %13.3f %13.3f %13.2f\n", nphot, mreg, x, y, z);
1810 (TVirtualMCApplication::Instance())->Stepping();
1811 }
3a625972 1812}
829fb838 1813