1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
18 Revision 1.13 2003/01/31 14:01:51 morsch
20 - Getters related to geometry.
21 - Communication with run manager (event steering)
23 Revision 1.12 2003/01/31 12:18:53 morsch
24 Corrected indices. (E. Futo)
26 Revision 1.9 2002/12/06 12:41:29 morsch
27 Mess from last merge cleaned up.
29 Revision 1.8 2002/12/06 12:28:44 morsch
30 Region to media mapping corrected and improved.
32 Revision 1.7 2002/12/06 12:21:32 morsch
33 User stepping methods added (E. Futo)
35 Revision 1.6 2002/11/21 18:40:06 iglez2
38 Revision 1.5 2002/11/07 17:59:10 iglez2
39 Included the geometry through geant4_vmc/FLUGG
41 Revision 1.4 2002/11/04 16:00:46 iglez2
42 The conversion between ID and PDG now uses Fluka routines and arrays which is more consistent.
44 Revision 1.3 2002/10/22 15:12:14 alibrary
45 Introducing Riostream.h
47 Revision 1.2 2002/10/14 14:57:40 hristov
48 Merging the VirtualMC branch to the main development branch (HEAD)
50 Revision 1.1.2.8 2002/10/08 16:33:17 iglez2
51 LSOUIT is set to true before the second call to flukam.
53 Revision 1.1.2.7 2002/10/08 09:30:37 iglez2
54 Solved stupid missing ;
56 Revision 1.1.2.6 2002/10/07 13:40:22 iglez2
57 First implementations of the PDG <--> Fluka Id conversion routines
59 Revision 1.1.2.5 2002/09/26 16:26:03 iglez2
61 Call to gAlice->Generator()->Generate()
63 Revision 1.1.2.4 2002/09/26 13:22:23 iglez2
64 Naive implementation of ProcessRun and ProcessEvent
65 Opening/Closing of input file (fInputFileName) with FORTRAN unit 5 before/after the first call to flukam inside Init()
67 Revision 1.1.2.3 2002/09/20 15:35:51 iglez2
68 Modification of LFDRTR. Value is passed to FLUKA !!!
70 Revision 1.1.2.2 2002/09/18 14:34:44 iglez2
71 Revised version with all pure virtual methods implemented
73 Revision 1.1.2.1 2002/07/24 08:49:41 alibrary
74 Adding TFluka to VirtualMC
76 Revision 1.1 2002/07/05 13:10:07 morsch
77 First commit of Fluka interface.
81 #include <Riostream.h>
83 #include "TClonesArray.h"
85 #include "TCallf77.h" //For the fortran calls
86 #include "Fdblprc.h" //(DBLPRC) fluka common
87 #include "Fepisor.h" //(EPISOR) fluka common
88 #include "Ffinuc.h" //(FINUC) fluka common
89 #include "Fiounit.h" //(IOUNIT) fluka common
90 #include "Fpaprop.h" //(PAPROP) fluka common
91 #include "Fpart.h" //(PART) fluka common
92 #include "Ftrackr.h" //(TRACKR) fluka common
93 #include "Fpaprop.h" //(PAPROP) fluka common
94 #include "Ffheavy.h" //(FHEAVY) fluka common
96 #include "TVirtualMC.h"
97 #include "TG4GeometryManager.h" //For the geometry management
98 #include "TG4DetConstruction.h" //For the detector construction
100 #include "FGeometryInit.hh"
101 #include "TLorentzVector.h"
102 #include "FlukaVolume.h"
104 // Fluka methods that may be needed.
106 # define flukam flukam_
107 # define fluka_openinp fluka_openinp_
108 # define fluka_closeinp fluka_closeinp_
109 # define mcihad mcihad_
110 # define mpdgha mpdgha_
112 # define flukam FLUKAM
113 # define fluka_openinp FLUKA_OPENINP
114 # define fluka_closeinp FLUKA_CLOSEINP
115 # define mcihad MCIHAD
116 # define mpdgha MPDGHA
122 // Prototypes for FLUKA functions
124 void type_of_call flukam(const int&);
125 void type_of_call fluka_openinp(const int&, DEFCHARA);
126 void type_of_call fluka_closeinp(const int&);
127 int type_of_call mcihad(const int&);
128 int type_of_call mpdgha(const int&);
132 // Class implementation for ROOT
137 //----------------------------------------------------------------------------
138 // TFluka constructors and destructors.
139 //____________________________________________________________________________
145 fCurrentFlukaRegion(-1)
148 // Default constructor
152 TFluka::TFluka(const char *title, Int_t verbosity)
153 :TVirtualMC("TFluka",title),
154 fVerbosityLevel(verbosity),
157 fCurrentFlukaRegion(-1)
159 if (fVerbosityLevel >=3)
160 cout << "==> TFluka::TFluka(" << title << ") constructor called." << endl;
163 // create geometry manager
164 if (fVerbosityLevel >=2)
165 cout << "\t* Creating G4 Geometry manager..." << endl;
166 fGeometryManager = new TG4GeometryManager();
167 if (fVerbosityLevel >=2)
168 cout << "\t* Creating G4 Detector..." << endl;
169 fDetector = new TG4DetConstruction();
170 FGeometryInit* geominit = FGeometryInit::GetInstance();
172 geominit->setDetConstruction(fDetector);
174 cerr << "ERROR: Could not create FGeometryInit!" << endl;
175 cerr << " Exiting!!!" << endl;
179 if (fVerbosityLevel >=3)
180 cout << "<== TFluka::TFluka(" << title << ") constructor called." << endl;
182 fVolumeMediaMap = new TClonesArray("FlukaVolume",1000);
188 if (fVerbosityLevel >=3)
189 cout << "==> TFluka::~TFluka() destructor called." << endl;
191 delete fGeometryManager;
192 fVolumeMediaMap->Delete();
193 delete fVolumeMediaMap;
196 if (fVerbosityLevel >=3)
197 cout << "<== TFluka::~TFluka() destructor called." << endl;
201 //_____________________________________________________________________________
202 // TFluka control methods
203 //____________________________________________________________________________
204 void TFluka::Init() {
205 if (fVerbosityLevel >=3)
206 cout << "==> TFluka::Init() called." << endl;
208 if (fVerbosityLevel >=2)
209 cout << "\t* Changing lfdrtr = (" << (GLOBAL.lfdrtr?'T':'F')
210 << ") in fluka..." << endl;
211 GLOBAL.lfdrtr = true;
213 if (fVerbosityLevel >=2)
214 cout << "\t* Opening file " << fInputFileName << endl;
215 const char* fname = fInputFileName;
216 fluka_openinp(lunin, PASSCHARA(fname));
218 if (fVerbosityLevel >=2)
219 cout << "\t* Calling flukam..." << endl;
222 if (fVerbosityLevel >=2)
223 cout << "\t* Closing file " << fInputFileName << endl;
224 fluka_closeinp(lunin);
226 if (fVerbosityLevel >=3)
227 cout << "<== TFluka::Init() called." << endl;
233 void TFluka::FinishGeometry() {
235 // Build-up table with region to medium correspondance
239 if (fVerbosityLevel >=3)
240 cout << "==> TFluka::FinishGeometry() called." << endl;
242 // fGeometryManager->Ggclos();
244 FGeometryInit* flugg = FGeometryInit::GetInstance();
246 fMediaByRegion = new Int_t[fNVolumes+2];
247 for (Int_t i = 0; i < fNVolumes; i++)
249 FlukaVolume* vol = dynamic_cast<FlukaVolume*>((*fVolumeMediaMap)[i]);
250 TString volName = vol->GetName();
251 Int_t media = vol->GetMedium();
252 printf("Finish Geometry: volName, media %d %s %d \n", i, volName.Data(), media);
253 strcpy(tmp, volName.Data());
255 flugg->SetMediumFromName(tmp, media, i+1);
256 fMediaByRegion[i] = media;
259 flugg->BuildMediaMap();
261 if (fVerbosityLevel >=3)
262 cout << "<== TFluka::FinishGeometry() called." << endl;
265 void TFluka::BuildPhysics() {
266 if (fVerbosityLevel >=3)
267 cout << "==> TFluka::BuildPhysics() called." << endl;
270 if (fVerbosityLevel >=3)
271 cout << "<== TFluka::BuildPhysics() called." << endl;
274 void TFluka::ProcessEvent() {
275 if (fVerbosityLevel >=3)
276 cout << "==> TFluka::ProcessEvent() called." << endl;
277 fApplication->GeneratePrimaries();
278 EPISOR.lsouit = true;
280 if (fVerbosityLevel >=3)
281 cout << "<== TFluka::ProcessEvent() called." << endl;
285 void TFluka::ProcessRun(Int_t nevent) {
286 if (fVerbosityLevel >=3)
287 cout << "==> TFluka::ProcessRun(" << nevent << ") called."
290 if (fVerbosityLevel >=2) {
291 cout << "\t* GLOBAL.fdrtr = " << (GLOBAL.lfdrtr?'T':'F') << endl;
292 cout << "\t* Calling flukam again..." << endl;
294 fApplication->InitGeometry();
295 fApplication->BeginEvent();
297 fApplication->FinishEvent();
298 if (fVerbosityLevel >=3)
299 cout << "<== TFluka::ProcessRun(" << nevent << ") called."
304 //_____________________________________________________________________________
305 // methods for building/management of geometry
306 //____________________________________________________________________________
307 // functions from GCONS
308 void TFluka::Gfmate(Int_t imat, char *name, Float_t &a, Float_t &z,
309 Float_t &dens, Float_t &radl, Float_t &absl,
310 Float_t* ubuf, Int_t& nbuf) {
312 fGeometryManager->Gfmate(imat, name, a, z, dens, radl, absl, ubuf, nbuf);
315 void TFluka::Gfmate(Int_t imat, char *name, Double_t &a, Double_t &z,
316 Double_t &dens, Double_t &radl, Double_t &absl,
317 Double_t* ubuf, Int_t& nbuf) {
319 fGeometryManager->Gfmate(imat, name, a, z, dens, radl, absl, ubuf, nbuf);
322 // detector composition
323 void TFluka::Material(Int_t& kmat, const char* name, Double_t a,
324 Double_t z, Double_t dens, Double_t radl, Double_t absl,
325 Float_t* buf, Int_t nwbuf) {
328 ->Material(kmat, name, a, z, dens, radl, absl, buf, nwbuf);
330 void TFluka::Material(Int_t& kmat, const char* name, Double_t a,
331 Double_t z, Double_t dens, Double_t radl, Double_t absl,
332 Double_t* buf, Int_t nwbuf) {
335 ->Material(kmat, name, a, z, dens, radl, absl, buf, nwbuf);
338 void TFluka::Mixture(Int_t& kmat, const char *name, Float_t *a,
339 Float_t *z, Double_t dens, Int_t nlmat, Float_t *wmat) {
342 ->Mixture(kmat, name, a, z, dens, nlmat, wmat);
344 void TFluka::Mixture(Int_t& kmat, const char *name, Double_t *a,
345 Double_t *z, Double_t dens, Int_t nlmat, Double_t *wmat) {
348 ->Mixture(kmat, name, a, z, dens, nlmat, wmat);
351 void TFluka::Medium(Int_t& kmed, const char *name, Int_t nmat,
352 Int_t isvol, Int_t ifield, Double_t fieldm, Double_t tmaxfd,
353 Double_t stemax, Double_t deemax, Double_t epsil,
354 Double_t stmin, Float_t* ubuf, Int_t nbuf) {
357 ->Medium(kmed, name, nmat, isvol, ifield, fieldm, tmaxfd, stemax, deemax,
358 epsil, stmin, ubuf, nbuf);
360 void TFluka::Medium(Int_t& kmed, const char *name, Int_t nmat,
361 Int_t isvol, Int_t ifield, Double_t fieldm, Double_t tmaxfd,
362 Double_t stemax, Double_t deemax, Double_t epsil,
363 Double_t stmin, Double_t* ubuf, Int_t nbuf) {
366 ->Medium(kmed, name, nmat, isvol, ifield, fieldm, tmaxfd, stemax, deemax,
367 epsil, stmin, ubuf, nbuf);
370 void TFluka::Matrix(Int_t& krot, Double_t thetaX, Double_t phiX,
371 Double_t thetaY, Double_t phiY, Double_t thetaZ,
375 ->Matrix(krot, thetaX, phiX, thetaY, phiY, thetaZ, phiZ);
378 void TFluka::Gstpar(Int_t itmed, const char *param, Double_t parval) {
380 fGeometryManager->Gstpar(itmed, param, parval);
383 // functions from GGEOM
384 Int_t TFluka::Gsvolu(const char *name, const char *shape, Int_t nmed,
385 Float_t *upar, Int_t np) {
387 // fVolumeMediaMap[TString(name)] = nmed;
388 printf("TFluka::Gsvolu() name = %s, nmed = %d\n", name, nmed);
390 TClonesArray &lvols = *fVolumeMediaMap;
391 new(lvols[fNVolumes++])
392 FlukaVolume(name, nmed);
393 return fGeometryManager->Gsvolu(name, shape, nmed, upar, np);
395 Int_t TFluka::Gsvolu(const char *name, const char *shape, Int_t nmed,
396 Double_t *upar, Int_t np) {
398 TClonesArray &lvols = *fVolumeMediaMap;
399 new(lvols[fNVolumes++])
400 FlukaVolume(name, nmed);
402 return fGeometryManager->Gsvolu(name, shape, nmed, upar, np);
405 void TFluka::Gsdvn(const char *name, const char *mother, Int_t ndiv,
408 // The medium of the daughter is the one of the mother
409 Int_t volid = TFluka::VolId(mother);
410 Int_t med = TFluka::VolId2Mate(volid);
411 TClonesArray &lvols = *fVolumeMediaMap;
412 new(lvols[fNVolumes++])
413 FlukaVolume(name, med);
414 fGeometryManager->Gsdvn(name, mother, ndiv, iaxis);
417 void TFluka::Gsdvn2(const char *name, const char *mother, Int_t ndiv,
418 Int_t iaxis, Double_t c0i, Int_t numed) {
420 TClonesArray &lvols = *fVolumeMediaMap;
421 new(lvols[fNVolumes++])
422 FlukaVolume(name, numed);
423 fGeometryManager->Gsdvn2(name, mother, ndiv, iaxis, c0i, numed);
426 void TFluka::Gsdvt(const char *name, const char *mother, Double_t step,
427 Int_t iaxis, Int_t numed, Int_t ndvmx) {
429 TClonesArray &lvols = *fVolumeMediaMap;
430 new(lvols[fNVolumes++])
431 FlukaVolume(name, numed);
432 fGeometryManager->Gsdvt(name, mother, step, iaxis, numed, ndvmx);
435 void TFluka::Gsdvt2(const char *name, const char *mother, Double_t step,
436 Int_t iaxis, Double_t c0, Int_t numed, Int_t ndvmx) {
438 TClonesArray &lvols = *fVolumeMediaMap;
439 new(lvols[fNVolumes++])
440 FlukaVolume(name, numed);
441 fGeometryManager->Gsdvt2(name, mother, step, iaxis, c0, numed, ndvmx);
444 void TFluka::Gsord(const char *name, Int_t iax) {
446 fGeometryManager->Gsord(name, iax);
449 void TFluka::Gspos(const char *name, Int_t nr, const char *mother,
450 Double_t x, Double_t y, Double_t z, Int_t irot,
453 fGeometryManager->Gspos(name, nr, mother, x, y, z, irot, konly);
456 void TFluka::Gsposp(const char *name, Int_t nr, const char *mother,
457 Double_t x, Double_t y, Double_t z, Int_t irot,
458 const char *konly, Float_t *upar, Int_t np) {
460 fGeometryManager->Gsposp(name, nr, mother, x, y, z, irot, konly, upar, np);
462 void TFluka::Gsposp(const char *name, Int_t nr, const char *mother,
463 Double_t x, Double_t y, Double_t z, Int_t irot,
464 const char *konly, Double_t *upar, Int_t np) {
466 fGeometryManager->Gsposp(name, nr, mother, x, y, z, irot, konly, upar, np);
469 void TFluka::Gsbool(const char* onlyVolName, const char* manyVolName) {
471 fGeometryManager->Gsbool(onlyVolName, manyVolName);
474 void TFluka::SetCerenkov(Int_t itmed, Int_t npckov, Float_t *ppckov,
475 Float_t *absco, Float_t *effic, Float_t *rindex) {
477 fGeometryManager->SetCerenkov(itmed, npckov, ppckov, absco, effic, rindex);
479 void TFluka::SetCerenkov(Int_t itmed, Int_t npckov, Double_t *ppckov,
480 Double_t *absco, Double_t *effic, Double_t *rindex) {
482 fGeometryManager->SetCerenkov(itmed, npckov, ppckov, absco, effic, rindex);
486 void TFluka::WriteEuclid(const char* fileName, const char* topVol,
487 Int_t number, Int_t nlevel) {
489 fGeometryManager->WriteEuclid(fileName, topVol, number, nlevel);
494 //_____________________________________________________________________________
495 // methods needed by the stepping
496 //____________________________________________________________________________
498 Int_t TFluka::GetMedium() const {
500 // Get the medium number for the current fluka region
502 FGeometryInit* flugg = FGeometryInit::GetInstance();
503 return flugg->GetMedium(fCurrentFlukaRegion);
508 //____________________________________________________________________________
509 // ID <--> PDG transformations
510 //_____________________________________________________________________________
511 Int_t TFluka::IdFromPDG(Int_t pdg) const
514 // Return Fluka code from PDG and pseudo ENDF code
516 // MCIHAD() goes from pdg to fluka internal.
517 Int_t intfluka = mcihad(pdg);
518 // KPTOIP array goes from internal to official
519 return GetFlukaKPTOIP(intfluka);
522 Int_t TFluka::PDGFromId(Int_t id) const
525 // Return PDG code and pseudo ENDF code from Fluka code
527 //IPTOKP array goes from official to internal
529 printf("PDGFromId: Error id = 0");
533 Int_t intfluka = GetFlukaIPTOKP(id);
535 printf("PDGFromId: Error intfluka = 0");
539 //MPKDHA() goes from internal to PDG
540 return mpdgha(intfluka);
543 //_____________________________________________________________________________
544 // methods for step management
545 //____________________________________________________________________________
549 void TFluka::SetMaxStep(Double_t)
551 // SetMaxStep is dummy procedure in TFluka !
552 cout << "SetMaxStep is dummy procedure in TFluka !" << endl;
555 void TFluka::SetMaxNStep(Int_t)
557 // SetMaxNStep is dummy procedure in TFluka !
558 cout << "SetMaxNStep is dummy procedure in TFluka !" << endl;
561 void TFluka::SetUserDecay(Int_t)
563 // SetUserDecay is dummy procedure in TFluka !
564 cout << "SetUserDecay is dummy procedure in TFluka !" << endl;
568 // dynamic properties
570 void TFluka::TrackPosition(TLorentzVector& position) const
572 // Return the current position in the master reference frame of the
573 // track being transported
574 // TRACKR.atrack = age of the particle
575 // TRACKR.xtrack = x-position of the last point
576 // TRACKR.ytrack = y-position of the last point
577 // TRACKR.ztrack = z-position of the last point
578 position.SetX(TRACKR.xtrack[TRACKR.ntrack]);
579 position.SetY(TRACKR.ytrack[TRACKR.ntrack]);
580 position.SetZ(TRACKR.ztrack[TRACKR.ntrack]);
581 position.SetT(TRACKR.atrack);
584 void TFluka::TrackMomentum(TLorentzVector& momentum) const
586 // Return the direction and the momentum (GeV/c) of the track
587 // currently being transported
588 // TRACKR.ptrack = momentum of the particle (not always defined, if
589 // < 0 must be obtained from etrack)
590 // TRACKR.cx,y,ztrck = direction cosines of the current particle
591 // TRACKR.etrack = total energy of the particle
592 // TRACKR.jtrack = identity number of the particle
593 // PAPROP.am[TRACKR.jtrack] = particle mass in gev
594 if (TRACKR.ptrack >= 0) {
595 momentum.SetPx(TRACKR.ptrack*TRACKR.cxtrck);
596 momentum.SetPy(TRACKR.ptrack*TRACKR.cytrck);
597 momentum.SetPz(TRACKR.ptrack*TRACKR.cztrck);
598 momentum.SetE(TRACKR.etrack);
602 Double_t p = sqrt(TRACKR.etrack*TRACKR.etrack - PAPROP.am[TRACKR.jtrack+6]*PAPROP.am[TRACKR.jtrack+6]);
603 momentum.SetPx(p*TRACKR.cxtrck);
604 momentum.SetPy(p*TRACKR.cytrck);
605 momentum.SetPz(p*TRACKR.cztrck);
606 momentum.SetE(TRACKR.etrack);
611 Double_t TFluka::TrackStep() const
613 // Return the length in centimeters of the current step
614 // TRACKR.ctrack = total curved path
615 return TRACKR.ctrack;
618 Double_t TFluka::TrackLength() const
621 // This is the sum of substeps !!!
622 // TRACKR.ctrack = total curved path of the current step
623 // Sum of the substeps is identical to TRACKR.ctrack if the is no mag. field
624 // The sum of all step length starting from the beginning of the track
625 // for the time being returns only the length in centimeters of the current step
627 for ( Int_t j=0;j<TRACKR.ntrack;j++) {
628 sum +=TRACKR.ttrack[j];
633 Double_t TFluka::TrackTime() const
635 // Return the current time of flight of the track being transported
636 // TRACKR.atrack = age of the particle
637 return TRACKR.atrack;
640 Double_t TFluka::Edep() const
643 // if TRACKR.ntrack = 0, TRACKR.mtrack = 0:
644 // -->local energy deposition (the value and the point are not recorded in TRACKR)
645 // but in the variable "rull" of the procedure "endraw.cxx"
646 // if TRACKR.ntrack > 0, TRACKR.mtrack = 0:
647 // -->no energy loss along the track
648 // if TRACKR.ntrack > 0, TRACKR.mtrack > 0:
649 // -->energy loss distributed along the track
650 // TRACKR.dtrack = energy deposition of the jth deposition even
651 if (TRACKR.ntrack == 0 && TRACKR.mtrack == 0)
655 for ( Int_t j=0;j<TRACKR.mtrack;j++) {
656 sum +=TRACKR.dtrack[j];
662 Int_t TFluka::TrackPid() const
664 // Return the id of the particle transported
665 // TRACKR.jtrack = identity number of the particle
666 return PDGFromId(TRACKR.jtrack);
669 Double_t TFluka::TrackCharge() const
671 // Return charge of the track currently transported
672 // PAPROP.ichrge = electric charge of the particle
673 // TRACKR.jtrack = identity number of the particle
674 return PAPROP.ichrge[TRACKR.jtrack+6];
677 Double_t TFluka::TrackMass() const
679 // PAPROP.am = particle mass in GeV
680 // TRACKR.jtrack = identity number of the particle
681 return PAPROP.am[TRACKR.jtrack+6];
684 Double_t TFluka::Etot() const
686 // TRACKR.etrack = total energy of the particle
687 return TRACKR.etrack;
693 Bool_t TFluka::IsNewTrack() const
696 // True if the track is not at the boundary of the current volume
697 // Not true in some cases in bxdraw - to be solved
701 Bool_t TFluka::IsTrackInside() const
703 // True if the track is not at the boundary of the current volume
704 // In Fluka a step is always inside one kind of material
705 // If the step would go behind the region of one material,
706 // it will be shortened to reach only the boundary.
707 // Therefore IsTrackInside() is always true.
708 // Not true in some cases in bxdraw - to be solved
712 Bool_t TFluka::IsTrackEntering() const
714 // True if this is the first step of the track in the current volume
715 // Boundary- (X) crossing
716 // Icode = 19: boundary crossing - call from Kaskad
717 // Icode = 29: boundary crossing - call from Emfsco
718 // Icode = 39: boundary crossing - call from Kasneu
719 // Icode = 49: boundary crossing - call from Kashea
720 // Icode = 59: boundary crossing - call from Kasoph
725 fIcode == 59) return 1;
729 Bool_t TFluka::IsTrackExiting() const
731 // True if this is the last step of the track in the current volume
732 // Boundary- (X) crossing
733 // Icode = 19: boundary crossing - call from Kaskad
734 // Icode = 29: boundary crossing - call from Emfsco
735 // Icode = 39: boundary crossing - call from Kasneu
736 // Icode = 49: boundary crossing - call from Kashea
737 // Icode = 59: boundary crossing - call from Kasoph
742 fIcode == 59) return 1;
746 Bool_t TFluka::IsTrackOut() const
748 // True if the track is out of the setup
750 // Icode = 14: escape - call from Kaskad
751 // Icode = 23: escape - call from Emfsco
752 // Icode = 32: escape - call from Kasneu
753 // Icode = 40: escape - call from Kashea
754 // Icode = 51: escape - call from Kasoph
759 fIcode == 51) return 1;
763 Bool_t TFluka::IsTrackDisappeared() const
765 // means all inelastic interactions and decays
766 // fIcode from usdraw
767 if (fIcode == 101 || // inelastic interaction
768 fIcode == 102 || // particle decay
769 fIcode == 214 || // in-flight annihilation
770 fIcode == 215 || // annihilation at rest
771 fIcode == 217 || // pair production
772 fIcode == 221) return 1;
776 Bool_t TFluka::IsTrackStop() const
778 // True if the track energy has fallen below the threshold
779 // means stopped by signal or below energy threshold
780 // Icode = 12: stopping particle - call from Kaskad
781 // Icode = 15: time kill - call from Kaskad
782 // Icode = 21: below threshold, iarg=1 - call from Emfsco
783 // Icode = 22: below threshold, iarg=2 - call from Emfsco
784 // Icode = 24: time kill - call from Emfsco
785 // Icode = 31: below threshold - call from Kasneu
786 // Icode = 33: time kill - call from Kasneu
787 // Icode = 41: time kill - call from Kashea
788 // Icode = 52: time kill - call from Kasoph
797 fIcode == 52) return 1;
801 Bool_t TFluka::IsTrackAlive() const
803 // means not disappeared or not out
804 if (IsTrackDisappeared() || IsTrackOut() ) return 0;
812 Int_t TFluka::NSecondaries() const
813 // Number of secondary particles generated in the current step
814 // FINUC.np = number of secondaries except light and heavy ions
815 // FHEAVY.npheav = number of secondaries for light and heavy secondary ions
817 return FINUC.np + FHEAVY.npheav;
820 void TFluka::GetSecondary(Int_t isec, Int_t& particleId,
821 TLorentzVector& position, TLorentzVector& momentum)
823 if (isec >= 0 && isec < FINUC.np) {
824 // more fine condition depending on icode
842 particleId = PDGFromId(FINUC.kpart[isec]);
843 position.SetX(fXsco);
844 position.SetY(fYsco);
845 position.SetZ(fZsco);
846 position.SetT(TRACKR.atrack);
847 // position.SetT(TRACKR.atrack+FINUC.agesec[isec]); //not yet implem.
848 momentum.SetPx(FINUC.plr[isec]*FINUC.cxr[isec]);
849 momentum.SetPy(FINUC.plr[isec]*FINUC.cyr[isec]);
850 momentum.SetPz(FINUC.plr[isec]*FINUC.czr[isec]);
851 momentum.SetE(FINUC.tki[isec] + PAPROP.am[FINUC.kpart[isec]+6]);
853 if (isec >= FINUC.np && isec < FINUC.np + FHEAVY.npheav) {
854 Int_t jsec = isec - FINUC.np;
855 particleId = FHEAVY.kheavy[jsec]; // this is Fluka id !!!
856 position.SetX(fXsco);
857 position.SetY(fYsco);
858 position.SetZ(fZsco);
859 position.SetT(TRACKR.atrack);
860 // position.SetT(TRACKR.atrack+FHEAVY.agheav[jsec]); //not yet implem.
861 momentum.SetPx(FHEAVY.pheavy[jsec]*FHEAVY.cxheav[jsec]);
862 momentum.SetPy(FHEAVY.pheavy[jsec]*FHEAVY.cyheav[jsec]);
863 momentum.SetPz(FHEAVY.pheavy[jsec]*FHEAVY.czheav[jsec]);
864 if (FHEAVY.tkheav[jsec] >= 3 && FHEAVY.tkheav[jsec] <= 6)
865 momentum.SetE(FHEAVY.tkheav[jsec] + PAPROP.am[jsec+6]);
866 else if (FHEAVY.tkheav[jsec] > 6)
867 momentum.SetE(FHEAVY.tkheav[jsec] + FHEAVY.amnhea[jsec]); // to be checked !!!
871 TMCProcess TFluka::ProdProcess(Int_t isec) const
872 // Name of the process that has produced the secondary particles
873 // in the current step
875 const TMCProcess kIpNoProc = kPNoProcess;
876 const TMCProcess kIpPDecay = kPDecay;
877 const TMCProcess kIpPPair = kPPair;
878 //const TMCProcess kIpPPairFromPhoton = kPPairFromPhoton;
879 //const TMCProcess kIpPPairFromVirtualPhoton = kPPairFromVirtualPhoton;
880 const TMCProcess kIpPCompton = kPCompton;
881 const TMCProcess kIpPPhotoelectric = kPPhotoelectric;
882 const TMCProcess kIpPBrem = kPBrem;
883 //const TMCProcess kIpPBremFromHeavy = kPBremFromHeavy;
884 //const TMCProcess kIpPBremFromElectronOrPositron = kPBremFromElectronOrPositron;
885 const TMCProcess kIpPDeltaRay = kPDeltaRay;
886 //const TMCProcess kIpPMoller = kPMoller;
887 //const TMCProcess kIpPBhabha = kPBhabha;
888 const TMCProcess kIpPAnnihilation = kPAnnihilation;
889 //const TMCProcess kIpPAnnihilInFlight = kPAnnihilInFlight;
890 //const TMCProcess kIpPAnnihilAtRest = kPAnnihilAtRest;
891 const TMCProcess kIpPHadronic = kPHadronic;
892 const TMCProcess kIpPMuonNuclear = kPMuonNuclear;
893 const TMCProcess kIpPPhotoFission = kPPhotoFission;
894 const TMCProcess kIpPRayleigh = kPRayleigh;
895 // const TMCProcess kIpPCerenkov = kPCerenkov;
896 // const TMCProcess kIpPSynchrotron = kPSynchrotron;
898 Int_t mugamma = TRACKR.jtrack == 7 || TRACKR.jtrack == 10 || TRACKR.jtrack == 11;
899 if (fIcode == 102) return kIpPDecay;
900 else if (fIcode == 104 || fIcode == 217) return kIpPPair;
901 //else if (fIcode == 104) return kIpPairFromPhoton;
902 //else if (fIcode == 217) return kIpPPairFromVirtualPhoton;
903 else if (fIcode == 219) return kIpPCompton;
904 else if (fIcode == 221) return kIpPPhotoelectric;
905 else if (fIcode == 105 || fIcode == 208) return kIpPBrem;
906 //else if (fIcode == 105) return kIpPBremFromHeavy;
907 //else if (fIcode == 208) return kPBremFromElectronOrPositron;
908 else if (fIcode == 103 || fIcode == 400) return kIpPDeltaRay;
909 else if (fIcode == 210 || fIcode == 212) return kIpPDeltaRay;
910 //else if (fIcode == 210) return kIpPMoller;
911 //else if (fIcode == 212) return kIpPBhabha;
912 else if (fIcode == 214 || fIcode == 215) return kIpPAnnihilation;
913 //else if (fIcode == 214) return kIpPAnnihilInFlight;
914 //else if (fIcode == 215) return kIpPAnnihilAtRest;
915 else if (fIcode == 101) return kIpPHadronic;
916 else if (fIcode == 101) {
917 if (!mugamma) return kIpPHadronic;
918 else if (TRACKR.jtrack == 7) return kIpPPhotoFission;
919 else return kIpPMuonNuclear;
921 else if (fIcode == 225) return kIpPRayleigh;
922 // Fluka codes 100, 300 and 400 still to be investigasted
923 else return kIpNoProc;
926 //Int_t StepProcesses(TArrayI &proc) const
927 // Return processes active in the current step
929 //ck = total energy of the particl ????????????????
933 Int_t TFluka::VolId2Mate(Int_t id) const
936 // Returns the material number for a given volume ID
938 printf("VolId2Mate %d %d\n", id, fMediaByRegion[id]);
939 return fMediaByRegion[id-1];
942 const char* TFluka::VolName(Int_t id) const
945 // Returns the volume name for a given volume ID
947 FlukaVolume* vol = dynamic_cast<FlukaVolume*>((*fVolumeMediaMap)[id-1]);
948 const char* name = vol->GetName();
949 printf("VolName %d %s \n", id, name);
953 Int_t TFluka::VolId(const Text_t* volName) const
956 // Converts from volume name to volume ID.
957 // Time consuming. (Only used during set-up)
958 // Could be replaced by hash-table
962 for (i = 0; i < fNVolumes; i++)
964 FlukaVolume* vol = dynamic_cast<FlukaVolume*>((*fVolumeMediaMap)[i]);
965 TString name = vol->GetName();
966 strcpy(tmp, name.Data());
968 if (!strcmp(tmp, volName)) break;
976 Int_t TFluka::CurrentVolID(Int_t& copyNo) const
979 // Return the logical id and copy number corresponding to the current fluka region
981 int ir = fCurrentFlukaRegion;
982 int id = (FGeometryInit::GetInstance())->CurrentVolID(ir, copyNo);
983 printf("CurrentVolID: %d %d %d \n", ir, id, copyNo);
988 Int_t TFluka::CurrentVolOffID(Int_t off, Int_t& copyNo) const
991 // Return the logical id and copy number of off'th mother
992 // corresponding to the current fluka region
995 return CurrentVolID(copyNo);
997 int ir = fCurrentFlukaRegion;
998 int id = (FGeometryInit::GetInstance())->CurrentVolOffID(ir, off, copyNo);
1000 printf("CurrentVolOffID: %d %d %d \n", ir, id, copyNo);
1002 printf("CurrentVolOffID: Warning Mother not found !!!\n");
1007 const char* TFluka::CurrentVolName() const
1010 // Return the current volume name
1013 Int_t id = TFluka::CurrentVolID(copy);
1014 const char* name = TFluka::VolName(id);
1015 printf("CurrentVolumeName: %d %s \n", fCurrentFlukaRegion, name);
1019 const char* TFluka::CurrentVolOffName(Int_t off) const
1022 // Return the volume name of the off'th mother of the current volume
1025 Int_t id = TFluka::CurrentVolOffID(off, copy);
1026 const char* name = TFluka::VolName(id);
1027 printf("CurrentVolumeOffName: %d %s \n", fCurrentFlukaRegion, name);
1031 Int_t TFluka::CurrentMaterial(Float_t &a, Float_t &z,
1032 Float_t &dens, Float_t &radl, Float_t &absl) const
1035 // Return the current medium number
1038 Int_t id = TFluka::CurrentVolID(copy);
1039 Int_t med = TFluka::VolId2Mate(id);
1040 printf("CurrentMaterial: %d %d \n", fCurrentFlukaRegion, med);
1044 void TFluka::Gmtod(Float_t* xm, Float_t* xd, Int_t iflag)
1046 // Transforms a position from the world reference frame
1047 // to the current volume reference frame.
1049 // Geant3 desription:
1050 // ==================
1051 // Computes coordinates XD (in DRS)
1052 // from known coordinates XM in MRS
1053 // The local reference system can be initialized by
1054 // - the tracking routines and GMTOD used in GUSTEP
1055 // - a call to GMEDIA(XM,NUMED)
1056 // - a call to GLVOLU(NLEVEL,NAMES,NUMBER,IER)
1057 // (inverse routine is GDTOM)
1059 // If IFLAG=1 convert coordinates
1060 // IFLAG=2 convert direction cosinus
1063 Double_t xmD[3], xdD[3];
1064 xmD[0] = xm[0]; xmD[1] = xm[1]; xmD[2] = xm[2];
1065 (FGeometryInit::GetInstance())->Gmtod(xmD, xdD, iflag);
1066 xd[0] = xdD[0]; xd[1] = xdD[1]; xd[2] = xdD[2];
1070 void TFluka::Gmtod(Double_t* xm, Double_t* xd, Int_t iflag)
1072 // Transforms a position from the world reference frame
1073 // to the current volume reference frame.
1075 // Geant3 desription:
1076 // ==================
1077 // Computes coordinates XD (in DRS)
1078 // from known coordinates XM in MRS
1079 // The local reference system can be initialized by
1080 // - the tracking routines and GMTOD used in GUSTEP
1081 // - a call to GMEDIA(XM,NUMED)
1082 // - a call to GLVOLU(NLEVEL,NAMES,NUMBER,IER)
1083 // (inverse routine is GDTOM)
1085 // If IFLAG=1 convert coordinates
1086 // IFLAG=2 convert direction cosinus
1089 Double_t xmD[3], xdD[3];
1090 xdD[0] = xd[0]; xdD[1] = xd[1]; xdD[2] = xd[2];
1091 (FGeometryInit::GetInstance())->Gdtom(xmD, xdD, iflag);
1092 xm[0] = xmD[0]; xm[1] = xmD[1]; xm[2] = xmD[2];
1095 void TFluka::Gdtom(Float_t* xd, Float_t* xm, Int_t iflag)
1097 // Transforms a position from the current volume reference frame
1098 // to the world reference frame.
1100 // Geant3 desription:
1101 // ==================
1102 // Computes coordinates XM (Master Reference System
1103 // knowing the coordinates XD (Detector Ref System)
1104 // The local reference system can be initialized by
1105 // - the tracking routines and GDTOM used in GUSTEP
1106 // - a call to GSCMED(NLEVEL,NAMES,NUMBER)
1107 // (inverse routine is GMTOD)
1109 // If IFLAG=1 convert coordinates
1110 // IFLAG=2 convert direction cosinus
1116 void TFluka::Gdtom(Double_t* xd, Double_t* xm, Int_t iflag)
1118 // Transforms a position from the current volume reference frame
1119 // to the world reference frame.
1121 // Geant3 desription:
1122 // ==================
1123 // Computes coordinates XM (Master Reference System
1124 // knowing the coordinates XD (Detector Ref System)
1125 // The local reference system can be initialized by
1126 // - the tracking routines and GDTOM used in GUSTEP
1127 // - a call to GSCMED(NLEVEL,NAMES,NUMBER)
1128 // (inverse routine is GMTOD)
1130 // If IFLAG=1 convert coordinates
1131 // IFLAG=2 convert direction cosinus
1135 (FGeometryInit::GetInstance())->Gdtom(xm, xd, iflag);
1138 // ===============================================================
1139 void TFluka::FutoTest()
1141 Int_t icode, mreg, newreg, particleId;
1143 Double_t rull, xsco, ysco, zsco;
1144 TLorentzVector position, momentum;
1147 cout << " icode=" << icode << endl;
1149 cout << "TLorentzVector positionX=" << position.X()
1150 << "positionY=" << position.Y()
1151 << "positionZ=" << position.Z()
1152 << "timeT=" << position.T() << endl;
1153 cout << "TLorentzVector momentumX=" << momentum.X()
1154 << "momentumY=" << momentum.Y()
1155 << "momentumZ=" << momentum.Z()
1156 << "energyE=" << momentum.E() << endl;
1157 cout << "TrackPid=" << TrackPid() << endl;
1161 else if (icode > 0 && icode <= 5) {
1164 // medium = GetMedium();
1165 cout << " icode=" << icode
1167 // << " medium=" << medium
1169 TrackPosition(position);
1170 TrackMomentum(momentum);
1171 cout << "TLorentzVector positionX=" << position.X()
1172 << "positionY=" << position.Y()
1173 << "positionZ=" << position.Z()
1174 << "timeT=" << position.T() << endl;
1175 cout << "TLorentzVector momentumX=" << momentum.X()
1176 << "momentumY=" << momentum.Y()
1177 << "momentumZ=" << momentum.Z()
1178 << "energyE=" << momentum.E() << endl;
1179 cout << "TrackStep=" << TrackStep() << endl;
1180 cout << "TrackLength=" << TrackLength() << endl;
1181 cout << "TrackTime=" << TrackTime() << endl;
1182 cout << "Edep=" << Edep() << endl;
1183 cout << "TrackPid=" << TrackPid() << endl;
1184 cout << "TrackCharge=" << TrackCharge() << endl;
1185 cout << "TrackMass=" << TrackMass() << endl;
1186 cout << "Etot=" << Etot() << endl;
1187 cout << "IsNewTrack=" << IsNewTrack() << endl;
1188 cout << "IsTrackInside=" << IsTrackInside() << endl;
1189 cout << "IsTrackEntering=" << IsTrackEntering() << endl;
1190 cout << "IsTrackExiting=" << IsTrackExiting() << endl;
1191 cout << "IsTrackOut=" << IsTrackOut() << endl;
1192 cout << "IsTrackDisappeared=" << IsTrackDisappeared() << endl;
1193 cout << "IsTrackAlive=" << IsTrackAlive() << endl;
1196 else if((icode >= 10 && icode <= 15) ||
1197 (icode >= 20 && icode <= 24) ||
1198 (icode >= 30 && icode <= 33) ||
1199 (icode >= 40 && icode <= 41) ||
1200 (icode >= 50 && icode <= 52)) {
1203 // medium = GetMedium();
1208 cout << " icode=" << icode
1210 // << " medium=" << medium
1214 << " zsco=" << zsco << endl;
1215 TrackPosition(position);
1216 TrackMomentum(momentum);
1217 cout << "Edep=" << Edep() << endl;
1218 cout << "Etot=" << Etot() << endl;
1219 cout << "TrackPid=" << TrackPid() << endl;
1220 cout << "TrackCharge=" << TrackCharge() << endl;
1221 cout << "TrackMass=" << TrackMass() << endl;
1222 cout << "IsTrackOut=" << IsTrackOut() << endl;
1223 cout << "IsTrackDisappeared=" << IsTrackDisappeared() << endl;
1224 cout << "IsTrackStop=" << IsTrackStop() << endl;
1225 cout << "IsTrackAlive=" << IsTrackAlive() << endl;
1228 else if((icode >= 100 && icode <= 105) ||
1232 (icode >= 214 && icode <= 215) ||
1241 // medium = GetMedium();
1245 cout << " icode=" << icode
1247 // << " medium=" << medium
1250 << " zsco=" << zsco << endl;
1251 cout << "TrackPid=" << TrackPid() << endl;
1252 cout << "NSecondaries=" << NSecondaries() << endl;
1253 for (Int_t isec=0; isec< NSecondaries(); isec++) {
1254 TFluka::GetSecondary(isec, particleId, position, momentum);
1255 cout << "TLorentzVector positionX=" << position.X()
1256 << "positionY=" << position.Y()
1257 << "positionZ=" << position.Z()
1258 << "timeT=" << position.T() << endl;
1259 cout << "TLorentzVector momentumX=" << momentum.X()
1260 << "momentumY=" << momentum.Y()
1261 << "momentumZ=" << momentum.Z()
1262 << "energyE=" << momentum.E() << endl;
1263 cout << "TrackPid=" << particleId << endl;
1268 else if((icode == 19) ||
1274 // medium = GetMedium();
1275 newreg = GetNewreg();
1279 cout << " icode=" << icode
1281 // << " medium=" << medium
1282 << " newreg=" << newreg
1285 << " zsco=" << zsco << endl;
1288 // ====================================================================
1293 } // end of FutoTest