From fa3d1cc75b38f4b4e3564dd43297764826fdb92c Mon Sep 17 00:00:00 2001 From: morsch Date: Fri, 6 Dec 2002 12:21:32 +0000 Subject: [PATCH] User stepping methods added (E. Futo) --- TFluka/Fdrawcalls.h | 104 ++++++++++ TFluka/Ffheavy.h | 89 ++++++++ TFluka/Ffinuc.h | 80 +++++++ TFluka/Fpaprop.h | 14 ++ TFluka/Ftrackr.h | 130 ++++++++++++ TFluka/TFluka.cxx | 483 ++++++++++++++++++++++++++++++++++++++++++- TFluka/TFluka.h | 110 ++++++---- TFluka/bxdraw.cxx | 24 +++ TFluka/eedraw.cxx | 18 ++ TFluka/endraw.cxx | 23 +++ TFluka/libTFluka.pkg | 8 +- TFluka/mgdraw.cxx | 20 ++ TFluka/sodraw.cxx | 18 ++ TFluka/usdraw.cxx | 23 +++ 14 files changed, 1096 insertions(+), 48 deletions(-) create mode 100644 TFluka/Fdrawcalls.h create mode 100644 TFluka/Ffheavy.h create mode 100644 TFluka/Ffinuc.h create mode 100644 TFluka/Ftrackr.h create mode 100644 TFluka/bxdraw.cxx create mode 100644 TFluka/eedraw.cxx create mode 100644 TFluka/endraw.cxx create mode 100644 TFluka/mgdraw.cxx create mode 100644 TFluka/sodraw.cxx create mode 100644 TFluka/usdraw.cxx diff --git a/TFluka/Fdrawcalls.h b/TFluka/Fdrawcalls.h new file mode 100644 index 00000000000..b7111bc9637 --- /dev/null +++ b/TFluka/Fdrawcalls.h @@ -0,0 +1,104 @@ +#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; + diff --git a/TFluka/Ffheavy.h b/TFluka/Ffheavy.h new file mode 100644 index 00000000000..b51e113dfec --- /dev/null +++ b/TFluka/Ffheavy.h @@ -0,0 +1,89 @@ +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); +} diff --git a/TFluka/Ffinuc.h b/TFluka/Ffinuc.h new file mode 100644 index 00000000000..6f9501595b5 --- /dev/null +++ b/TFluka/Ffinuc.h @@ -0,0 +1,80 @@ +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); +} diff --git a/TFluka/Fpaprop.h b/TFluka/Fpaprop.h index 5a6786447e8..27c5db6ef1e 100644 --- a/TFluka/Fpaprop.h +++ b/TFluka/Fpaprop.h @@ -86,4 +86,18 @@ typedef struct { 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 diff --git a/TFluka/Ftrackr.h b/TFluka/Ftrackr.h new file mode 100644 index 00000000000..4b27e514c11 --- /dev/null +++ b/TFluka/Ftrackr.h @@ -0,0 +1,130 @@ +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 = 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];}; +} diff --git a/TFluka/TFluka.cxx b/TFluka/TFluka.cxx index 2d6cd10ee2f..b439dc5742f 100644 --- a/TFluka/TFluka.cxx +++ b/TFluka/TFluka.cxx @@ -15,6 +15,9 @@ /* $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 @@ -63,15 +66,20 @@ First commit of Fluka interface. #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 @@ -190,6 +198,7 @@ void TFluka::Init() { if (fVerbosityLevel >=3) cout << "<== TFluka::Init() called." << endl; + } void TFluka::FinishGeometry() { @@ -447,3 +456,473 @@ Int_t TFluka::PDGFromId(Int_t id) const 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= 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 + diff --git a/TFluka/TFluka.h b/TFluka/TFluka.h index ce72dd9608c..9bbaf47549d 100644 --- a/TFluka/TFluka.h +++ b/TFluka/TFluka.h @@ -191,52 +191,31 @@ class TFluka : public TVirtualMC { // 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 @@ -295,23 +274,65 @@ class TFluka : public TVirtualMC { // - 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 @@ -320,7 +341,10 @@ class TFluka : public TVirtualMC { std::vector fMediaByRegion; + ClassDef(TFluka,1) //C++ interface to Fluka montecarlo + + }; #endif //TFLUKA diff --git a/TFluka/bxdraw.cxx b/TFluka/bxdraw.cxx new file mode 100644 index 00000000000..90832c0e550 --- /dev/null +++ b/TFluka/bxdraw.cxx @@ -0,0 +1,24 @@ +#include +#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" + diff --git a/TFluka/eedraw.cxx b/TFluka/eedraw.cxx new file mode 100644 index 00000000000..04290f4073b --- /dev/null +++ b/TFluka/eedraw.cxx @@ -0,0 +1,18 @@ +#include +#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" + diff --git a/TFluka/endraw.cxx b/TFluka/endraw.cxx new file mode 100644 index 00000000000..0213740711e --- /dev/null +++ b/TFluka/endraw.cxx @@ -0,0 +1,23 @@ +#include +#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" + diff --git a/TFluka/libTFluka.pkg b/TFluka/libTFluka.pkg index d7caee3f8a8..d05bfaba667 100644 --- a/TFluka/libTFluka.pkg +++ b/TFluka/libTFluka.pkg @@ -1,9 +1,11 @@ # 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 diff --git a/TFluka/mgdraw.cxx b/TFluka/mgdraw.cxx new file mode 100644 index 00000000000..5d1d4dc172d --- /dev/null +++ b/TFluka/mgdraw.cxx @@ -0,0 +1,20 @@ +#include +#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" + diff --git a/TFluka/sodraw.cxx b/TFluka/sodraw.cxx new file mode 100644 index 00000000000..5ae563c2fe8 --- /dev/null +++ b/TFluka/sodraw.cxx @@ -0,0 +1,18 @@ +#include +#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" + diff --git a/TFluka/usdraw.cxx b/TFluka/usdraw.cxx new file mode 100644 index 00000000000..da0e61424de --- /dev/null +++ b/TFluka/usdraw.cxx @@ -0,0 +1,23 @@ +#include +#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" + -- 2.43.0