--- /dev/null
+#include "Rtypes.h"
+//
+// Structure saving parameters of Fluka procedures
+// mgdraw(Int_t& icode, Int_t& mreg)
+// bxdraw(Int_t& icode, Int_t& mreg, Int_t& newreg,
+// Double_t& xsco, Double_t& ysco, Double_t& zsco)
+// eedraw(Int_t& icode)
+// endraw(Int_t& icode, Int_t& mreg,
+// Double_t& rull, Double_t& xsco, Double_t& ysco, Double_t& zsco)
+// sogdraw
+// usdraw(Int_t& icode, Int_t& mreg,
+// Double_t& xsco, Double_t& ysco, Double_t& zsco)
+//
+//*---------------------- icode values ----------------------------------*
+//* *
+//* Icode = 1: call from Kaskad *
+//* Icode = 2: call from Emfsco *
+//* Icode = 3: call from Kasneu *
+//* Icode = 4: call from Kashea *
+//* Icode = 5: call from Kasoph *
+//* *
+//* Boundary- (X) crossing *
+//* ====================== *
+//* Icode = 1x: call from Kaskad *
+//* Icode = 19: boundary crossing *
+//* Icode = 2x: call from Emfsco *
+//* Icode = 29: boundary crossing *
+//* Icode = 3x: call from Kasneu *
+//* Icode = 39: boundary crossing *
+//* Icode = 4x: call from Kashea *
+//* Icode = 49: boundary crossing *
+//* Icode = 59: call from Kasoph *
+//* Icode = 59: boundary crossing *
+//* *
+//* ENergy deposition DRAWing: *
+//* ========================= *
+//* *
+//* Icode = 1x: call from Kaskad *
+//* 10: elastic interaction recoil *
+//* 11: inelastic interaction recoil *
+//* 12: stopping particle *
+//* 13: pseudo-neutron deposition *
+//* 14: escape *
+//* 15: time kill *
+//* Icode = 2x: call from Emfsco *
+//* 20: local energy deposition (i.e. photoelectric) *
+//* 21: below threshold, iarg=1 *
+//* 22: below threshold, iarg=2 *
+//* 23: escape *
+//* 24: time kill *
+//* Icode = 3x: call from Kasneu *
+//* 30: target recoil *
+//* 31: below threshold *
+//* 32: escape *
+//* 33: time kill *
+//* Icode = 4x: call from Kashea *
+//* 40: escape *
+//* 41: time kill *
+//* Icode = 5x: call from Kasoph *
+//* 50: optical photon absorption *
+//* 51: escape *
+//* 52: time kill *
+//* *
+//* USer dependent DRAWing: *
+//* ====================== *
+//* *
+//* Icode = 10x: call from Kaskad *
+//* 100: elastic interaction secondaries *
+//* 101: inelastic interaction secondaries *
+//* 102: particle decay secondaries *
+//* 103: delta ray generation secondaries *
+//* 104: pair production secondaries *
+//* 105: bremsstrahlung secondaries *
+//* Icode = 20x: call from Emfsco *
+//* 208: bremsstrahlung secondaries *
+//* 210: Moller secondaries *
+//* 212: Bhabha secondaries *
+//* 214: in-flight annihilation secondaries *
+//* 215: annihilation at rest secondaries *
+//* 217: pair production secondaries *
+//* 219: Compton scattering secondaries *
+//* 221: photoelectric secondaries *
+//* 225: Rayleigh scattering secondaries *
+//* Icode = 30x: call from Kasneu *
+//* 300: interaction secondaries *
+//* Icode = 40x: call from Kashea *
+//* 400: delta ray generation secondaries *
+//* For all interactions secondaries are put on FINUC common (kp=1,np)
+//* but for KASHEA delta ray generation where only the secondary *
+//* electron is present and stacked on STACK common for kp=lstack *
+//* *
+//*---------------------- end of icode values ---------------------------*
+//
+
+typedef struct {
+ Int_t icode;
+ Int_t mreg;
+ Int_t newreg;
+ Double_t rull;
+ Double_t xsco;
+ Double_t ysco;
+ Double_t zsco;
+} Fdrawcalls;
+
--- /dev/null
+extern "C" {
+//*$ create fheavy.add
+//*copy fheavy
+//*
+//*=== fheavy ===========================================================*
+//*
+//*----------------------------------------------------------------------*
+//* *
+//* include file: fheavy *
+//* *
+//* created on 5 april 1990 by alfredo ferrari, infn milan *
+//* *
+//* last change on 26-jul-97 by alfredo ferrari, infn milan *
+//* *
+//* included in the following subroutines or functions: not updated *
+//* *
+//* description of the common block(s) and variable(s) *
+//* *
+//* /fheavy/ is the storage for heavy secondaries created in the *
+//* nuclear evaporation *
+//* npheav = number of secondaries *
+//* kheavy(ip) = type of the secondary ip *
+//* ( 3 = deuteron, 4 = 3-h, 5 = 3-he, 6 = 4-he, *
+//* 7-12 = "heavy" fragment specified by ibheav and *
+//* icheav ) *
+//* cxheav(ip) = direction cosine of the secondary ip *
+//* with respect to x-axis *
+//* cyheav(ip) = direction cosine of the secondary ip *
+//* with respect to y-axis *
+//* czheav(ip) = direction cosine of the secondary ip *
+//* with respect to z-axis *
+//* tkheav(ip) = kinetic energy of secondary ip *
+//* pheavy(ip) = momentum of the secondary ip *
+//* wheavy(ip) = weight of the secondary ip *
+//* agheav(ip) = "age" of the secondary ip with respect to the *
+//* interaction time *
+//* amheav(kp) = atomic masses of the twelve types of evaporated *
+//* or fragmented or fissioned particles *
+//* amnhea(kp) = nuclear masses of the twelve types of evaporated *
+//* or fragmented or fissioned particles *
+//* bhheav(jp,kp) = (nuclear) binding energy of the jp_th hyperon of *
+//* the kp-type heavy particle *
+//* anheav(kp) = name of the kp-type heavy particle *
+//* icheav(kp) = charge of the kp-type heavy particle *
+//* ibheav(kp) = mass number of the kp-type heavy particle *
+//* imheav(kp) = isomeric state of the kp-type heavy particle *
+//* ihheav(kp) = number of hyperons of the kp-type heavy particle *
+//* khheav(jp,kp) = id of the jp_th hyperon of the kp-type heavy *
+//* particle *
+//* !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! *
+//* !!! there is now the possibility to produce up to 6 "heavy" !!!! *
+//* !!! fragments besides the residual nucleus recorded in !!!! *
+//* !!! resnuc: they are identified by indeces 7-12, of course !!!! *
+//* !!! the corresponding physical properties (z,a,m..) must be !!!! *
+//* !!! updated every time they are produced !!!! *
+//* !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! *
+//*----------------------------------------------------------------------*
+//*
+const Int_t mxheav = 100;
+const Int_t kxheav = 12;
+
+typedef struct {
+ Double_t cxheav[mxheav];
+ Double_t cyheav[mxheav];
+ Double_t czheav[mxheav];
+ Double_t tkheav[mxheav];
+ Double_t pheavy[mxheav];
+ Double_t wheavy[mxheav];
+ Double_t agheav[mxheav];
+ Double_t bhheav[kxheav][ihypmx];
+ Double_t amheav[kxheav];
+ Double_t amnhea[kxheav];
+ Int_t kheavy[mxheav];
+ Int_t icheav[kxheav];
+ Int_t ibheav[kxheav];
+ Int_t imheav[kxheav];
+ Int_t ihheav[kxheav];
+ Int_t khheav[kxheav][ihypmx];
+ Int_t npheav;
+} fheavyCommon;
+#define FHEAVY COMMON_BLOCK(FHEAVY,fheavy)
+COMMON_BLOCK_DEF(fheavyCommon,FHEAVY);
+
+typedef struct {
+ Char_t anheav[kxheav][8];
+} fheavcCommon;
+#define FHEAVC COMMON_BLOCK(FHEAVC,fheavc)
+COMMON_BLOCK_DEF(fheavcCommon,FHEAVC);
+}
--- /dev/null
+extern "C" {
+//*$ create finuc.add
+//*copy finuc
+//*
+//*=== finuc ============================================================*
+//*
+//*----------------------------------------------------------------------*
+//* *
+//* include file: finuc (new version of old finuc of fluka86) *
+//* *
+//* !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! *
+//* !!!! s e e a l s o i n c l u d e f i l e !!!! *
+//* !!!! f i n u c 2 !!!! *
+//* !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! *
+//* *
+//* created on 20 january 1996 by alfredo ferrari & paola sala *
+//* infn - milan *
+//* *
+//* last change on 26-jul-97 by alfredo ferrari *
+//* *
+//* included in the following subroutines or functions: not updated *
+//* *
+//* description of the common block(s) and variable(s) *
+//* *
+//* /finuc/ is the storage for secondaries created in event *
+//* np = number of secondaries *
+//* kpart (ip) = type of the secondary ip *
+//* cxr (ip) = direction cosine of the secondary ip *
+//* with respect to x-axis *
+//* cyr (ip) = direction cosine of the secondary ip *
+//* with respect to y-axis *
+//* czr (ip) = direction cosine of the secondary ip *
+//* with respect to z-axis *
+//* cxrpol (ip) = direction cosine of the secondary ip polarization *
+//* with respect to x-axis *
+//* cyrpol (ip) = direction cosine of the secondary ip polarization *
+//* with respect to y-axis *
+//* czrpol (ip) = direction cosine of the secondary ip polarization *
+//* with respect to z-axis *
+//* tki (ip) = kinetic energy of secondary ip *
+//* plr (ip) = momentum of the secondary ip *
+//* wei (ip) = weight of the secondary ip *
+//* agesec (ip) = "age" of the secondary ip with respect to the *
+//* interaction time *
+//* tv = excitation energy *
+//* tvcms = actual excitation energy of the residual nucleus *
+//* tvrecl = recoil kinetic energy of the residual nucleus *
+//* tvheav = recoil kinetic energies of heavy (2-h, 3-h, 3-he, *
+//* 4-he) fragments after evaporation *
+//* tvbind = approximate energy wasted in nuclear binding *
+//* effects (not yet operational) *
+//* *
+//*----------------------------------------------------------------------*
+//*
+const Int_t mxp = 999;
+//*
+
+typedef struct {
+ Double_t cxr[mxp];
+ Double_t cyr[mxp];
+ Double_t czr[mxp];
+ Double_t cxrpol[mxp];
+ Double_t cyrpol[mxp];
+ Double_t czrpol[mxp];
+ Double_t tki[mxp];
+ Double_t plr[mxp];
+ Double_t wei[mxp];
+ Double_t agesec[mxp];
+ Double_t tv;
+ Double_t tvcms;
+ Double_t tvrecl;
+ Double_t tvheav;
+ Double_t tvbind;
+ Int_t np0;
+ Int_t np;
+ Int_t kpart[mxp];
+} finucCommon;
+#define FINUC COMMON_BLOCK(FINUC,finuc)
+COMMON_BLOCK_DEF(finucCommon,FINUC);
+}
COMMON_BLOCK_DEF(chpprpCommon,CHPPRP);
}
+//Get functions
+//inline Double_t GetFlukaAM(unsigned int i) {return PAPROP.am[i+6];}
+//inline Double_t GetFlukaAMDISC(unsigned int i) {return PAPROP.amdisc[i+6];}
+//inline Double_t GetFlukaTHALH(unsigned int i) {return PAPROP.thalf[i+6];}
+//inline Double_t GetFlukaBIASDC(unsigned int i) {return PAPROP.biasdc[i+6];}
+//inline Double_t GetFlukaBIASIN(unsigned int i) {return PAPROP.biasin[i+6];}
+//inline Int_t GetFlukaICHRGE(unsigned int i) {return PAPROP.ichrge[i+6];}
+//inline Int_t GetFlukaIBARCH(unsigned int i) {return PAPROP.ibarch[i+6];}
+//inline Int_t GetFlukaIJDISC(unsigned int i) {return PAPROP.ijdisc[i+6];}
+//inline Int_t GetFlukaJSPINP(unsigned int i) {return PAPROP.jspinp[i+6];}
+//inline Int_t GetFlukaIPARTY(unsigned int i) {return PAPROP.iparty[i+6];}
+//inline Int_t GetFlukaIPARID(unsigned int i) {return PAPROP.iparid[i+6];}
+//inline Int_t GetFlukaLHADRO(unsigned int i) {return PAPROP.lhadro[i+6];}
+//inline Int_t GetFlukaLBSDCY(unsigned int i) {return PAPROP.lbsdcy[i+6];}
#endif
--- /dev/null
+extern "C" {
+//*$ create trackr.add
+//*copy trackr
+//* *
+//*=== trackr ===========================================================*
+//* *
+//*----------------------------------------------------------------------*
+//* *
+//* tracks recording by alfredo ferrari, infn - milan *
+//* *
+//* last change 31 january 2001 by alfredo ferrari *
+//* *
+//* included in : *
+//* electr *
+//* emfsco *
+//* kaskad (new version) *
+//* kashea *
+//* kasneu *
+//* geoden (new version) *
+//* mageas *
+//* magmov *
+//* magnew *
+//* move *
+//* photon *
+//* usrsco *
+//* *
+//* ntrack = number of track segments *
+//* mtrack = number of energy deposition events along the track *
+//* 0 < i < ntrack *
+//* xtrack = end x-point of the ith track segment *
+//* ytrack = end y-point of the ith track segment *
+//* ztrack = end z-point of the ith track segment *
+//* 1 < i < ntrack *
+//* ttrack = length of the ith track segment *
+//* 1 < j < mtrack *
+//* dtrack = energy deposition of the jth deposition event *
+//* *
+//* jtrack = identity number of the particle *
+//* etrack = total energy of the particle *
+//* ptrack = momentum of the particle (not always defined, if *
+//* < 0 must be obtained from etrack) *
+//* cx,y,ztrck = direction cosines of the current particle *
+//* cx,y,ztrpl = polarization cosines of the current particle *
+//* wtrack = weight of the particle *
+//* wscrng = scoring weight: it can differ from wtrack if some *
+//* biasing techniques are used (for example inelastic *
+//* interaction length biasing) *
+//* ctrack = total curved path *
+//* zfftrk = <z_eff> of the particle *
+//* zfrttk = actual z_eff of the particle *
+//* atrack = age of the particle *
+//* akshrt = kshrt amplitude for k0/k0bar *
+//* aklong = klong amplitude for k0/k0bar *
+//* wninou = neutron algebraic balance of interactions (both *
+//* for "high" energy particles and "low" energy *
+//* neutrons) *
+//* spausr = user defined spare variables for the current *
+//* particle *
+//* sttrck = macroscopic total cross section for low energy *
+//* neutron collisions *
+//* satrck = macroscopic absorption cross section for low energy*
+//* neutron collisions (it can be negative for pnab>1) *
+//* ktrack = if > 0 neutron group of the particle (neutron) *
+//* *
+//* ntrack > 0, mtrack > 0 : energy loss distributed along the *
+//* track *
+//* ntrack > 0, mtrack = 0 : no energy loss along the track *
+//* ntrack = 0, mtrack = 0 : local energy deposition (the *
+//* value and the point are not re- *
+//* corded in trackr) *
+//* mmtrck = flag recording the material index for low energy *
+//* neutron collisions *
+//* lt1trk = initial lattice cell of the current track *
+//* (or lattice cell for a point energy deposition) *
+//* lt2trk = final lattice cell of the current track *
+//* ihspnt = current geometry history pointer (not set if -1) *
+//* ltrack = flag recording the generation number *
+//* llouse = user defined flag for the current particle *
+//* ispusr = user defined spare flags for the current particle *
+//* lfsssc = logical flag for inelastic interactions ending with*
+//* fission (used also for low energy neutrons) *
+//* *
+//*----------------------------------------------------------------------*
+//* *
+const Int_t mxtrck = 2500;
+
+typedef struct {
+ Double_t xtrack[mxtrck+1];
+ Double_t ytrack[mxtrck+1];
+ Double_t ztrack[mxtrck+1];
+ Double_t ttrack[mxtrck];
+ Double_t dtrack[mxtrck];
+ Double_t etrack;
+ Double_t ptrack;
+ Double_t cxtrck;
+ Double_t cytrck;
+ Double_t cztrck;
+ Double_t wtrack;
+ Double_t cxtrpl;
+ Double_t cytrpl;
+ Double_t cztrpl;
+ Double_t zfftrk;
+ Double_t zfrttk;
+ Double_t atrack;
+ Double_t ctrack;
+ Double_t akshrt;
+ Double_t aklong;
+ Double_t wscrng;
+ Double_t wninou;
+ Double_t spausr[mkbmx1];
+ Double_t sttrck;
+ Double_t satrck;
+ Int_t ntrack;
+ Int_t mtrack;
+ Int_t jtrack;
+ Int_t ktrack;
+ Int_t mmtrck;
+ Int_t lt1trk;
+ Int_t lt2trk;
+ Int_t ihspnt;
+ Int_t ltrack;
+ Int_t llouse;
+ Int_t ispusr[mkbmx2];
+ Int_t lfsssc;
+} trackrCommon;
+#define TRACKR COMMON_BLOCK(TRACKR,trackr)
+COMMON_BLOCK_DEF(trackrCommon,TRACKR);
+//static union { Double_t spause; Double_t spausr[0];};
+//static union { Int_t ispuse; Int_t ispusr[0];};
+}
/*
$Log$
+Revision 1.6 2002/11/21 18:40:06 iglez2
+Media handling added
+
Revision 1.5 2002/11/07 17:59:10 iglez2
Included the geometry through geant4_vmc/FLUGG
#include "TFluka.h"
#include "TCallf77.h" //For the fortran calls
#include "Fdblprc.h" //(DBLPRC) fluka common
-#include "Fiounit.h" //(IOUNIT) fluka common
#include "Fepisor.h" //(EPISOR) fluka common
+#include "Ffinuc.h" //(FINUC) fluka common
+#include "Fiounit.h" //(IOUNIT) fluka common
+#include "Fpaprop.h" //(PAPROP) fluka common
#include "Fpart.h" //(PART) fluka common
-#include "TVirtualMC.h"
+#include "Ftrackr.h" //(TRACKR) fluka common
+#include "Ffheavy.h" //(FHEAVY) fluka common
+#include "TVirtualMC.h"
#include "TG4GeometryManager.h" //For the geometry management
#include "TG4DetConstruction.h" //For the detector construction
#include "FGeometryInit.hh"
+#include "TLorentzVector.h"
// Fluka methods that may be needed.
#ifndef WIN32
if (fVerbosityLevel >=3)
cout << "<== TFluka::Init() called." << endl;
+
}
void TFluka::FinishGeometry() {
return mpdgha(intfluka);
}
+
+
+
+//_____________________________________________________________________________
+// methods for step management
+//____________________________________________________________________________
+//
+// dynamic properties
+//
+void TFluka::TrackPosition(TLorentzVector& position) const
+{
+// Return the current position in the master reference frame of the
+// track being transported
+// TRACKR.atrack = age of the particle
+// TRACKR.xtrack = x-position of the last point
+// TRACKR.ytrack = y-position of the last point
+// TRACKR.ztrack = z-position of the last point
+ position.SetX(TRACKR.xtrack[TRACKR.ntrack]);
+ position.SetY(TRACKR.ytrack[TRACKR.ntrack]);
+ position.SetZ(TRACKR.ztrack[TRACKR.ntrack]);
+ position.SetT(TRACKR.atrack);
+}
+
+void TFluka::TrackMomentum(TLorentzVector& momentum) const
+{
+// Return the direction and the momentum (GeV/c) of the track
+// currently being transported
+// TRACKR.ptrack = momentum of the particle (not always defined, if
+// < 0 must be obtained from etrack)
+// TRACKR.cx,y,ztrck = direction cosines of the current particle
+// TRACKR.etrack = total energy of the particle
+// TRACKR.jtrack = identity number of the particle
+// PAPROP.am[TRACKR.jtrack] = particle mass in gev
+ if (TRACKR.ptrack >= 0) {
+ momentum.SetPx(TRACKR.ptrack*TRACKR.cxtrck);
+ momentum.SetPy(TRACKR.ptrack*TRACKR.cytrck);
+ momentum.SetPz(TRACKR.ptrack*TRACKR.cztrck);
+ momentum.SetE(TRACKR.etrack);
+ return;
+ }
+ else {
+ Double_t p = sqrt(TRACKR.etrack*TRACKR.etrack - PAPROP.am[TRACKR.jtrack]*PAPROP.am[TRACKR.jtrack]);
+ momentum.SetPx(p*TRACKR.cxtrck);
+ momentum.SetPy(p*TRACKR.cytrck);
+ momentum.SetPz(p*TRACKR.cztrck);
+ momentum.SetE(TRACKR.etrack);
+ return;
+ }
+}
+
+Double_t TFluka::TrackStep() const
+{
+// Return the length in centimeters of the current step
+// TRACKR.ctrack = total curved path
+ return TRACKR.ctrack;
+}
+
+Double_t TFluka::TrackLength() const
+{
+// It is wrong
+// should be the sum of all steps starting from the beginning of the track
+// for the time being returns only the length in centimeters of the current step
+ return TRACKR.ctrack;
+}
+
+Double_t TFluka::TrackTime() const
+{
+// Return the current time of flight of the track being transported
+// TRACKR.atrack = age of the particle
+ return TRACKR.atrack;
+}
+
+Double_t TFluka::Edep() const
+{
+// Energy deposition
+// if TRACKR.ntrack = 0, TRACKR.mtrack = 0:
+// -->local energy deposition (the value and the point are not recorded in TRACKR)
+// but in the variable "rull" of the procedure "endraw.cxx"
+// if TRACKR.ntrack > 0, TRACKR.mtrack = 0:
+// -->no energy loss along the track
+// if TRACKR.ntrack > 0, TRACKR.mtrack > 0:
+// -->energy loss distributed along the track
+// TRACKR.dtrack = energy deposition of the jth deposition even
+ if (TRACKR.ntrack == 0 && TRACKR.mtrack == 0)
+ return fRull;
+ else {
+ Double_t sum = 0;
+ for ( Int_t j=0;j<TRACKR.mtrack;j++) {
+ sum +=TRACKR.dtrack[j];
+ }
+ return sum;
+ }
+}
+
+Int_t TFluka::TrackPid() const
+{
+// Return the id of the particle transported
+// TRACKR.jtrack = identity number of the particle
+ return PDGFromId(TRACKR.jtrack);
+}
+
+Double_t TFluka::TrackCharge() const
+{
+// Return charge of the track currently transported
+// PAPROP.ichrge = electric charge of the particle
+ return PAPROP.ichrge[TRACKR.jtrack+6];
+}
+
+Double_t TFluka::TrackMass() const
+{
+// PAPROP.am = particle mass in GeV
+ return PAPROP.am[TRACKR.jtrack+6];
+}
+
+Double_t TFluka::Etot() const
+{
+// TRACKR.etrack = total energy of the particle
+ return TRACKR.etrack;
+}
+
+//
+// track status
+//
+Bool_t TFluka::IsNewTrack() const
+{
+// ???????????????,
+// True if the track is not at the boundary of the current volume
+// Not true in some cases in bxdraw - to be solved
+ return 1;
+}
+
+Bool_t TFluka::IsTrackInside() const
+{
+// True if the track is not at the boundary of the current volume
+// In Fluka a step is always inside one kind of material
+// If the step would go behind the region of one material,
+// it will be shortened to reach only the boundary.
+// Therefore IsTrackInside() is always true.
+// Not true in some cases in bxdraw - to be solved
+ return 1;
+}
+
+Bool_t TFluka::IsTrackEntering() const
+{
+// True if this is the first step of the track in the current volume
+// Boundary- (X) crossing
+// Icode = 19: boundary crossing - call from Kaskad
+// Icode = 29: boundary crossing - call from Emfsco
+// Icode = 39: boundary crossing - call from Kasneu
+// Icode = 49: boundary crossing - call from Kashea
+// Icode = 59: boundary crossing - call from Kasoph
+ if (fIcode == 19 ||
+ fIcode == 29 ||
+ fIcode == 39 ||
+ fIcode == 49 ||
+ fIcode == 59) return 1;
+ else return 0;
+}
+
+Bool_t TFluka::IsTrackExiting() const
+{
+// True if this is the last step of the track in the current volume
+// Boundary- (X) crossing
+// Icode = 19: boundary crossing - call from Kaskad
+// Icode = 29: boundary crossing - call from Emfsco
+// Icode = 39: boundary crossing - call from Kasneu
+// Icode = 49: boundary crossing - call from Kashea
+// Icode = 59: boundary crossing - call from Kasoph
+ if (fIcode == 19 ||
+ fIcode == 29 ||
+ fIcode == 39 ||
+ fIcode == 49 ||
+ fIcode == 59) return 1;
+ else return 0;
+}
+
+Bool_t TFluka::IsTrackOut() const
+{
+// True if the track is out of the setup
+// means escape
+// Icode = 14: escape - call from Kaskad
+// Icode = 23: escape - call from Emfsco
+// Icode = 32: escape - call from Kasneu
+// Icode = 40: escape - call from Kashea
+// Icode = 51: escape - call from Kasoph
+ if (fIcode == 14 ||
+ fIcode == 23 ||
+ fIcode == 32 ||
+ fIcode == 40 ||
+ fIcode == 51) return 1;
+ else return 0;
+}
+
+Bool_t TFluka::IsTrackDisappeared() const
+{
+// means all inelastic interactions and decays
+// fIcode from usdraw
+ if (fIcode == 101 || // inelastic interaction
+ fIcode == 102 || // particle decay
+ fIcode == 214 || // in-flight annihilation
+ fIcode == 215 || // annihilation at rest
+ fIcode == 217 || // pair production
+ fIcode == 221) return 1;
+ else return 0;
+}
+
+Bool_t TFluka::IsTrackStop() const
+{
+// True if the track energy has fallen below the threshold
+// means stopped by signal or below energy threshold
+// Icode = 12: stopping particle - call from Kaskad
+// Icode = 15: time kill - call from Kaskad
+// Icode = 21: below threshold, iarg=1 - call from Emfsco
+// Icode = 22: below threshold, iarg=2 - call from Emfsco
+// Icode = 24: time kill - call from Emfsco
+// Icode = 31: below threshold - call from Kasneu
+// Icode = 33: time kill - call from Kasneu
+// Icode = 41: time kill - call from Kashea
+// Icode = 52: time kill - call from Kasoph
+ if (fIcode == 12 ||
+ fIcode == 15 ||
+ fIcode == 21 ||
+ fIcode == 22 ||
+ fIcode == 24 ||
+ fIcode == 31 ||
+ fIcode == 33 ||
+ fIcode == 41 ||
+ fIcode == 52) return 1;
+ else return 0;
+}
+
+Bool_t TFluka::IsTrackAlive() const
+{
+// means not disappeared or not out
+ if (IsTrackDisappeared() || IsTrackOut() ) return 0;
+ else return 1;
+}
+
+//
+// secondaries
+//
+
+Int_t TFluka::NSecondaries() const
+// Number of secondary particles generated in the current step
+// fIcode = 100 = elastic interaction
+// fIcode = 101 = inelastic interaction
+// fIcode = 102 = particle decay
+// fIcode = 103 = delta ray
+// fIcode = 104 = pair production
+// fIcode = 105 = bremsstrahlung
+{
+ if (fIcode >= 100 && fIcode <= 105)
+ return FINUC.np + FHEAVY.npheav;
+ else
+ return -1;
+}
+
+void TFluka::GetSecondary(Int_t isec, Int_t& particleId,
+ TLorentzVector& position, TLorentzVector& momentum)
+//
+// fIcode = 100 = elastic interaction
+// fIcode = 101 = inelastic interaction
+// fIcode = 102 = particle decay
+// fIcode = 103 = delta ray
+// fIcode = 104 = pair production
+// fIcode = 105 = bremsstrahlung
+{
+ if (fIcode >= 100 && fIcode <= 105) {
+ if (isec >= 0 && isec < FINUC.np) {
+ particleId = PDGFromId(FINUC.kpart[isec]);
+ position.SetX(fXsco);
+ position.SetY(fYsco);
+ position.SetZ(fZsco);
+ position.SetT(FINUC.agesec[isec]);
+ momentum.SetPx(FINUC.plr[isec]*FINUC.cxr[isec]);
+ momentum.SetPy(FINUC.plr[isec]*FINUC.cyr[isec]);
+ momentum.SetPz(FINUC.plr[isec]*FINUC.czr[isec]);
+ momentum.SetE(FINUC.tki[isec] + PAPROP.am[FINUC.kpart[isec]+6]);
+ }
+ if (isec >= FINUC.np && isec < FINUC.np + FHEAVY.npheav) {
+ Int_t jsec = isec - FINUC.np;
+ particleId = FHEAVY.kheavy[jsec]; // this is Fluka id !!!
+ position.SetX(fXsco);
+ position.SetY(fYsco);
+ position.SetZ(fZsco);
+ position.SetT(FHEAVY.agheav[jsec]);
+ momentum.SetPx(FHEAVY.pheavy[jsec]*FHEAVY.cxheav[jsec]);
+ momentum.SetPy(FHEAVY.pheavy[jsec]*FHEAVY.cyheav[jsec]);
+ momentum.SetPz(FHEAVY.pheavy[jsec]*FHEAVY.czheav[jsec]);
+ if (FHEAVY.tkheav[jsec] >= 3 && FHEAVY.tkheav[jsec] <= 6)
+ momentum.SetE(FHEAVY.tkheav[jsec] + PAPROP.am[jsec+6]);
+ else if (FHEAVY.tkheav[jsec] > 6)
+ momentum.SetE(FHEAVY.tkheav[jsec] + FHEAVY.amnhea[jsec]); // to be checked !!!
+ }
+ }
+}
+
+//TMCProcess ProdProcess(Int_t isec) const
+// Name of the process that has produced the secondary particles
+// in the current step
+//{
+// will come from FINUC when called from USDRAW
+//}
+
+//Int_t StepProcesses(TArrayI &proc) const
+// Return processes active in the current step
+//{
+//ck = total energy of the particl ????????????????
+//}
+
+
+// ===============================================================
+void TFluka::FutoTest()
+{
+ Int_t icode, mreg, newreg, particleId;
+// Int_t medium;
+ Double_t rull, xsco, ysco, zsco;
+ TLorentzVector position, momentum;
+ icode = GetIcode();
+ if (icode == 0) {
+ cout << " icode=" << icode << endl;
+ /*
+ cout << "TLorentzVector positionX=" << position.X()
+ << "positionY=" << position.Y()
+ << "positionZ=" << position.Z()
+ << "timeT=" << position.T() << endl;
+ cout << "TLorentzVector momentumX=" << momentum.X()
+ << "momentumY=" << momentum.Y()
+ << "momentumZ=" << momentum.Z()
+ << "energyE=" << momentum.E() << endl;
+ cout << "TrackPid=" << TrackPid() << endl;
+ */
+ }
+
+ else if (icode > 0 && icode <= 5) {
+// mgdraw
+ mreg = GetMreg();
+// medium = GetMedium();
+ cout << " icode=" << icode
+ << " mreg=" << mreg
+// << " medium=" << medium
+ << endl;
+ TrackPosition(position);
+ TrackMomentum(momentum);
+ cout << "TLorentzVector positionX=" << position.X()
+ << "positionY=" << position.Y()
+ << "positionZ=" << position.Z()
+ << "timeT=" << position.T() << endl;
+ cout << "TLorentzVector momentumX=" << momentum.X()
+ << "momentumY=" << momentum.Y()
+ << "momentumZ=" << momentum.Z()
+ << "energyE=" << momentum.E() << endl;
+ cout << "TrackStep=" << TrackStep() << endl;
+ cout << "TrackLength=" << TrackLength() << endl;
+ cout << "TrackTime=" << TrackTime() << endl;
+ cout << "Edep=" << Edep() << endl;
+ cout << "TrackPid=" << TrackPid() << endl;
+ cout << "TrackCharge=" << TrackCharge() << endl;
+ cout << "TrackMass=" << TrackMass() << endl;
+ cout << "Etot=" << Etot() << endl;
+ cout << "IsNewTrack=" << IsNewTrack() << endl;
+ cout << "IsTrackInside=" << IsTrackInside() << endl;
+ cout << "IsTrackEntering=" << IsTrackEntering() << endl;
+ cout << "IsTrackExiting=" << IsTrackExiting() << endl;
+ cout << "IsTrackOut=" << IsTrackOut() << endl;
+ cout << "IsTrackDisappeared=" << IsTrackDisappeared() << endl;
+ cout << "IsTrackAlive=" << IsTrackAlive() << endl;
+ }
+
+ else if((icode >= 10 && icode <= 15) ||
+ (icode >= 20 && icode <= 24) ||
+ (icode >= 30 && icode <= 33) ||
+ (icode >= 40 && icode <= 41) ||
+ (icode >= 50 && icode <= 52)) {
+// endraw
+ mreg = GetMreg();
+// medium = GetMedium();
+ rull = GetRull();
+ xsco = GetXsco();
+ ysco = GetYsco();
+ zsco = GetZsco();
+ cout << " icode=" << icode
+ << " mreg=" << mreg
+// << " medium=" << medium
+ << " rull=" << rull
+ << " xsco=" << xsco
+ << " ysco=" << ysco
+ << " zsco=" << zsco << endl;
+ TrackPosition(position);
+ TrackMomentum(momentum);
+ cout << "Edep=" << Edep() << endl;
+ cout << "Etot=" << Etot() << endl;
+ cout << "TrackPid=" << TrackPid() << endl;
+ cout << "TrackCharge=" << TrackCharge() << endl;
+ cout << "TrackMass=" << TrackMass() << endl;
+ cout << "IsTrackOut=" << IsTrackOut() << endl;
+ cout << "IsTrackDisappeared=" << IsTrackDisappeared() << endl;
+ cout << "IsTrackStop=" << IsTrackStop() << endl;
+ cout << "IsTrackAlive=" << IsTrackAlive() << endl;
+ }
+
+ else if((icode >= 100 && icode <= 105) ||
+ (icode == 208) ||
+ (icode == 210) ||
+ (icode == 212) ||
+ (icode >= 214 && icode <= 215) ||
+ (icode == 217) ||
+ (icode == 219) ||
+ (icode == 221) ||
+ (icode == 225) ||
+ (icode == 300) ||
+ (icode == 400)) {
+// usdraw
+ mreg = GetMreg();
+// medium = GetMedium();
+ xsco = GetXsco();
+ ysco = GetYsco();
+ zsco = GetZsco();
+ cout << " icode=" << icode
+ << " mreg=" << mreg
+// << " medium=" << medium
+ << " xsco=" << xsco
+ << " ysco=" << ysco
+ << " zsco=" << zsco << endl;
+ cout << "TrackPid=" << TrackPid() << endl;
+ cout << "NSecondaries=" << NSecondaries() << endl;
+ for (Int_t isec=0; isec< NSecondaries(); isec++) {
+//void TFluka::GetSecondary(Int_t isec, Int_t& particleId,
+// TLorentzVector& position, TLorentzVector& momentum)
+ TFluka::GetSecondary(isec, particleId, position, momentum);
+ cout << "TLorentzVector positionX=" << position.X()
+ << "positionY=" << position.Y()
+ << "positionZ=" << position.Z()
+ << "timeT=" << position.T() << endl;
+ cout << "TLorentzVector momentumX=" << momentum.X()
+ << "momentumY=" << momentum.Y()
+ << "momentumZ=" << momentum.Z()
+ << "energyE=" << momentum.E() << endl;
+ cout << "TrackPid=" << particleId << endl;
+
+ }
+ }
+
+ else if((icode == 19) ||
+ (icode == 29) ||
+ (icode == 39) ||
+ (icode == 49) ||
+ (icode == 59)) {
+ mreg = GetMreg();
+// medium = GetMedium();
+ newreg = GetNewreg();
+ xsco = GetXsco();
+ ysco = GetYsco();
+ zsco = GetZsco();
+ cout << " icode=" << icode
+ << " mreg=" << mreg
+// << " medium=" << medium
+ << " newreg=" << newreg
+ << " xsco=" << xsco
+ << " ysco=" << ysco
+ << " zsco=" << zsco << endl;
+ }
+//
+// ====================================================================
+//
+
+
+
+} // end of FutoTest
+
// tracking particle
// dynamic properties
- virtual void TrackPosition(TLorentzVector& position) const
- {printf("WARNING: TrackPosition not yet implemented !\n");}
- virtual void TrackMomentum(TLorentzVector& momentum) const
- {printf("WARNING: TrackMomentum not yet implemented !\n");}
- virtual Double_t TrackStep() const
- {printf("WARNING: TrackStep not yet implemented !\n"); return -1.;}
- virtual Double_t TrackLength() const
- {printf("WARNING: TrackLength not yet implemented !\n"); return -1.;}
- virtual Double_t TrackTime() const
- {printf("WARNING: TrackTime not yet implemented !\n"); return -1.;}
- virtual Double_t Edep() const
- {printf("WARNING: Edep not yet implemented !\n"); return -1.;}
+ virtual void TrackPosition(TLorentzVector& position) const;
+ virtual void TrackMomentum(TLorentzVector& momentum) const;
+ virtual Double_t TrackStep() const;
+ virtual Double_t TrackLength() const;
+ virtual Double_t TrackTime() const;
+ virtual Double_t Edep() const;
// static properties
- virtual Int_t TrackPid() const
- {printf("WARNING: TrackPid not yet implemented !\n"); return -1;}
- virtual Double_t TrackCharge() const
- {printf("WARNING: TrackCharge not yet implemented !\n"); return -1.;}
- virtual Double_t TrackMass() const
- {printf("WARNING: TrackMass not yet implemented !\n"); return -1.;}
- virtual Double_t Etot() const
- {printf("WARNING: Etot not yet implemented !\n"); return -1.;}
-
+ virtual Int_t TrackPid() const;
+ virtual Double_t TrackCharge() const;
+ virtual Double_t TrackMass() const;
+ virtual Double_t Etot() const;
// track status
- virtual Bool_t IsNewTrack() const
- {printf("WARNING: IsNewTrack not yet implemented !\n"); return 0;}
- virtual Bool_t IsTrackInside() const
- {printf("WARNING: IsTrackInside not yet implemented !\n"); return 0;}
- virtual Bool_t IsTrackEntering() const
- {printf("WARNING: IsTrackEntering not yet implemented !\n"); return 0;}
- virtual Bool_t IsTrackExiting() const
- {printf("WARNING: IsTrackExiting not yet implemented !\n"); return 0;}
- virtual Bool_t IsTrackOut() const
- {printf("WARNING: IsTrackOut not yet implemented !\n"); return 0;}
- virtual Bool_t IsTrackDisappeared() const
- {printf("WARNING: IsTrackDisappeared not yet implemented !\n"); return 0;}
- virtual Bool_t IsTrackStop() const
- {printf("WARNING: IsTrackStop not yet implemented !\n"); return 0;}
- virtual Bool_t IsTrackAlive() const
- {printf("WARNING: IsTrackAlive not yet implemented !\n"); return 0;}
+ virtual Bool_t IsNewTrack() const;
+ virtual Bool_t IsTrackInside() const;
+ virtual Bool_t IsTrackEntering() const;
+ virtual Bool_t IsTrackExiting() const;
+ virtual Bool_t IsTrackOut() const;
+ virtual Bool_t IsTrackDisappeared() const;
+ virtual Bool_t IsTrackStop() const;
+ virtual Bool_t IsTrackAlive() const;
// secondaries
- virtual Int_t NSecondaries() const
- {printf("WARNING: NSecondaries not yet implemented !\n"); return -1;}
+ virtual Int_t NSecondaries() const ;
virtual void GetSecondary(Int_t isec, Int_t& particleId,
- TLorentzVector& position, TLorentzVector& momentum)
- {printf("WARNING: GetSecondary not yet implemented !\n");}
+ TLorentzVector& position, TLorentzVector& momentum);
virtual TMCProcess ProdProcess(Int_t isec) const
{printf("WARNING: StepProcesses not yet implemented !\n"); return kPNoProcess;}
virtual Int_t StepProcesses(TArrayI &proc) const
// - Verbosity level
Int_t GetVerbosityLevel() const {return fVerbosityLevel;}
void SetVerbosityLevel(Int_t l) {fVerbosityLevel = l;}
+
+ // - Fluka Draw procedures formal parameters
+ Int_t GetIcode() const {return fIcode;}
+ void SetIcode(Int_t l) {fIcode = l;}
+
+ Int_t GetMreg() const {return fCurrentFlukaRegion;}
+ void SetMreg(Int_t l) {fCurrentFlukaRegion = l;}
+
+ Int_t GetNewreg() const {return fNewreg;}
+ void SetNewreg(Int_t l) {fNewreg = l;}
+
+ Double_t GetRull() const {return fRull;}
+ void SetRull(Double_t r) {fRull = r;}
+
+ Double_t GetXsco() const {return fXsco;}
+ void SetXsco(Double_t x) {fXsco = x;}
+
+ Double_t GetYsco() const {return fYsco;}
+ void SetYsco(Double_t y) {fYsco = y;}
+
+ Double_t GetZsco() const {return fZsco;}
+ void SetZsco(Double_t z) {fZsco = z;}
+
void SetCurrentFlukaRegion(Int_t reg) {fCurrentFlukaRegion=reg;}
Int_t GetCurrentFlukaRegion() const {return fCurrentFlukaRegion;}
+ //
+ // test
+ // ------------------------------------------------
+ //
+ virtual void FutoTest() ;
+
private:
TFluka(const TFluka &mc){}
TFluka & operator=(const TFluka &) {return (*this);}
protected:
Int_t fVerbosityLevel; //Verbosity level (0 lowest - 3 highest)
- TString fInputFileName; //Name of the input file (f.e. alice.inp)
-
+
+ TString fInputFileName; //Name of the input file (f.e. mu.inp)
+
+ Int_t fIcode; //Fluka Draw procedures formal parameter
+// Int_t fMreg; //Fluka Draw procedures formal parameter
+ Int_t fNewreg; //Fluka Draw procedures formal parameter
+ Double_t fRull; //Fluka Draw procedures formal parameter
+ Double_t fXsco; //Fluka Draw procedures formal parameter
+ Double_t fYsco; //Fluka Draw procedures formal parameter
+ Double_t fZsco; //Fluka Draw procedures formal parameter
+
//Geometry through Geant4 for the time being!!!
TG4GeometryManager* fGeometryManager; //geometry manager
TG4DetConstruction* fDetector; //Detector
+
+
+
+
//Index of fluka region at each step
Int_t fCurrentFlukaRegion;
//Map between volume name and media indices
std::vector<Int_t> fMediaByRegion;
+
ClassDef(TFluka,1) //C++ interface to Fluka montecarlo
+
+
};
#endif //TFLUKA
--- /dev/null
+#include <Riostream.h>
+#include "AliRun.h"
+#include "TFluka.h"
+#ifndef WIN32
+# define bxdraw bxdraw_
+#else
+# define bxdraw BXDRAW
+#endif
+extern "C" {
+void bxdraw(Int_t& icode, Int_t& mreg, Int_t& newreg,
+ Double_t& xsco, Double_t& ysco, Double_t& zsco)
+{
+ ((TFluka*) gMC)->SetIcode(icode);
+ ((TFluka*) gMC)->SetMreg(mreg);
+ ((TFluka*) gMC)->SetNewreg(newreg);
+ ((TFluka*) gMC)->SetXsco(xsco);
+ ((TFluka*) gMC)->SetYsco(ysco);
+ ((TFluka*) gMC)->SetZsco(zsco);
+ cout << endl << " !!! I am in bxdraw - calling gAlice->Stepping()" << endl;
+ ((TFluka*) gMC)->FutoTest();
+// gAlice->Stepping();
+} // end of bxdraw
+} // end of extern "C"
+
--- /dev/null
+#include <Riostream.h>
+#include "AliRun.h"
+#include "TFluka.h"
+#ifndef WIN32
+# define eedraw eedraw_
+#else
+# define eedraw EEDRAW
+#endif
+extern "C" {
+void eedraw(Int_t& icode)
+{
+ ((TFluka*) gMC)->SetIcode(icode);
+ cout << endl << " !!! I am in eedraw - calling gAlice->Stepping()" << endl;
+ ((TFluka*) gMC)->FutoTest();
+// gAlice->Stepping();
+} // end of eedraw
+} // end of extern "C"
+
--- /dev/null
+#include <Riostream.h>
+#include "AliRun.h"
+#include "TFluka.h"
+#ifndef WIN32
+# define endraw endraw_
+#else
+# define endraw ENDRAW
+#endif
+extern "C" {
+void endraw(Int_t& icode, Int_t& mreg, Double_t& rull, Double_t& xsco, Double_t& ysco, Double_t& zsco)
+{
+ ((TFluka*) gMC)->SetIcode(icode);
+ ((TFluka*) gMC)->SetMreg(mreg);
+ ((TFluka*) gMC)->SetRull(rull);
+ ((TFluka*) gMC)->SetXsco(xsco);
+ ((TFluka*) gMC)->SetYsco(ysco);
+ ((TFluka*) gMC)->SetZsco(zsco);
+ cout << endl << " !!! I am in endraw - calling gAlice->Stepping()" << endl;
+ ((TFluka*) gMC)->FutoTest();
+// gAlice->Stepping();
+} // end of endraw
+} // end of extern "C"
+
# Sources
-SRCS:= TFluka.cxx source.cxx
-FSRCS:= FLUKA_input.f
+SRCS:= TFluka.cxx source.cxx mgdraw.cxx bxdraw.cxx eedraw.cxx endraw.cxx sodraw.cxx usdraw.cxx
+#FSRCS:= FLUKA_input.f mgdraw.f
+FSRCS:= FLUKA_input.f
# Headers
-HDRS:= TFluka.h
+
+HDRS:= TFluka.h
# ROOT Dictionary
DHDR:= TFlukaLinkDef.h
--- /dev/null
+#include <Riostream.h>
+#include "AliRun.h"
+#include "TFluka.h"
+#ifndef WIN32
+# define mgdraw mgdraw_
+#else
+# define mgdraw MGDRAW
+#endif
+
+extern "C" {
+void mgdraw(Int_t& icode, Int_t& mreg)
+{
+ ((TFluka*) gMC)->SetIcode(icode);
+ ((TFluka*) gMC)->SetMreg(mreg);
+ cout << endl << " !!! I am in mgdraw - calling gAlice->Stepping()" << endl;
+ ((TFluka*) gMC)->FutoTest();
+// gAlice->Stepping();
+} // end of mgdraw
+} // end of extern "C"
+
--- /dev/null
+#include <Riostream.h>
+#include "AliRun.h"
+#include "TFluka.h"
+#ifndef WIN32
+# define sodraw sodraw_
+#else
+# define sodraw SODRAW
+#endif
+extern "C" {
+void sodraw()
+{
+ ((TFluka*) gMC)->SetIcode(0);
+ cout << endl << " !!! I am in sodraw - calling gAlice->Stepping()" << endl;
+ ((TFluka*) gMC)->FutoTest();
+// gAlice->Stepping();
+} // end of sodraw
+} // end of extern "C"
+
--- /dev/null
+#include <Riostream.h>
+#include "AliRun.h"
+#include "TFluka.h"
+#ifndef WIN32
+# define usdraw usdraw_
+#else
+# define usdraw USDRAW
+#endif
+extern "C" {
+void usdraw(Int_t& icode, Int_t& mreg,
+ Double_t& xsco, Double_t& ysco, Double_t& zsco)
+{
+ ((TFluka*) gMC)->SetIcode(icode);
+ ((TFluka*) gMC)->SetMreg(mreg);
+ ((TFluka*) gMC)->SetXsco(xsco);
+ ((TFluka*) gMC)->SetYsco(ysco);
+ ((TFluka*) gMC)->SetZsco(zsco);
+ cout << endl << " !!! I am in usdraw - calling gAlice->Stepping()" << endl;
+ ((TFluka*) gMC)->FutoTest();
+// gAlice->Stepping();
+} // end of usdraw
+} // end of extern "C"
+
git://git.uio.no / u / mrichter / AliRoot.git / commitdiff