From 5929ad29610764495bfcdd12cf6e3b09be95de10 Mon Sep 17 00:00:00 2001 From: morsch Date: Mon, 26 Jan 2004 16:29:58 +0000 Subject: [PATCH] Updates needed for geant4.6 and fluka2004. --- TFluka/Fbeam.h | 54 ++++++++---------- TFluka/Fdimpar.h | 8 ++- TFluka/Femfstk.h | 1 + TFluka/Ffheavy.h | 2 +- TFluka/Ffinuc.h | 10 +++- TFluka/Fpaprop.h | 3 +- TFluka/Fstack.h | 3 + TFluka/Ftrackr.h | 4 ++ TFluka/TFluka.cxx | 130 +++++++++++++++++++++++++++++-------------- TFluka/bxdraw.cxx | 2 +- TFluka/source_bg.cxx | 8 ++- 11 files changed, 142 insertions(+), 83 deletions(-) diff --git a/TFluka/Fbeam.h b/TFluka/Fbeam.h index 5617bb739a1..fe94ae14416 100644 --- a/TFluka/Fbeam.h +++ b/TFluka/Fbeam.h @@ -41,33 +41,29 @@ extern "C" { //* tinpz = direction cosine of the beam polariz. with respect to* //* z-axis * //* polfra = polarization fraction * -//* nforce = number of the region of forced interaction * -//* xfor = x-coord. of the starting point of the region nforce * -//* yfor = y-coord. of the starting point of the region nforce * -//* zfor = z-coord. of the starting point of the region nforce * -//* disfor = thickness of the region nforce in cm * -//* wfor = relative weight of the particle due to forcing * +//* beawei = weight of the beam particles * +//* bmaxis(j,i) = j_th component of the i_th axis used to define the* +//* conventional x,y,z beam reference frame * +//*!!!!! ATTENTION: in C++ it is the component bmaxis(i,j) !!!!! * //* ijbeam = beam particle type (see btype in /paprop/) * //* ijhion = heavy ion type if ijbeam = -2 * -//* ipbite = flag describing the shape of the momentum * -//* distribution of the beam * -//* 0=rectangular, 1=gaussian * -//* idiv = flag describing the shape of the angular * -//* divergence distribution of the beam * -//* 0=rectangular, 1=gaussian * -//* ixspot = flag describing the shape of the spatial * -//* distribution of the beam spot in x-direction * -//* 0=rectangular, 1=gaussian * -//* iyspot = flag describing the shape of the spatial * -//* distribution of the beam spot in y-direction * -//* 0=rectangular, 1=gaussian * -//* beawei = weight of the beam particles * +//* ldpgss = true for a gaussian momentum distribution of the * +//* beam particles, false for a rectangular one * +//* ldvgss = true for a gaussian angular divergence distribution * +//* of the beam particles, false for a rectangular one * +//* ldxgss = true for a gaussian spatial distribution of the beam * +//* spot in the x-direction, false for a rectangular one * +//* ldygss = true for a gaussian spatial distribution of the beam * +//* spot in the y-direction, false for a rectangular one * //* lbeamc = flag for an annular beam * //* lpperp = flag for polar. perp. to the beam direction * -//* lpfrac = flag for interpreting the polar. fraction * +//* lpfrac = flag for interpreting the polar. fraction +//* lbaxis = logical flag for using a beam axis frame different * +//* from the standard one * //* * //*----------------------------------------------------------------------* -// LOGICAL LBEAMC, LPPERP, LPFRAC +// LBEAMC, LPPERP, LPFRAC, LDPGSS, LDVGSS, LDXGSS, LDYGSS, +// & LBAXIS typedef struct { Double_t pbeam; @@ -86,21 +82,17 @@ typedef struct { Double_t tinpz; Double_t polfra; Double_t beawei; - Double_t xfor; - Double_t yfor; - Double_t zfor; - Double_t disfor; - Double_t wfor; + Double_t bmaxis[3][3]; Int_t ijbeam; Int_t ijhion; - Int_t ipbite; - Int_t idiv; - Int_t ixspot; - Int_t iyspot; - Int_t nforce; + Int_t ldpgss; + Int_t ldvgss; + Int_t ldxgss; + Int_t ldygss; Int_t lbeamc; Int_t lpperp; Int_t lpfrac; + Int_t lbaxis; } beamCommon; #define BEAM COMMON_BLOCK(BEAM,beam) COMMON_BLOCK_DEF(beamCommon,BEAM); diff --git a/TFluka/Fdimpar.h b/TFluka/Fdimpar.h index f7283806c4f..98dffaf7380 100644 --- a/TFluka/Fdimpar.h +++ b/TFluka/Fdimpar.h @@ -1,5 +1,5 @@ #ifndef FDIMPAR_H -#define FDIMPAR_H 1 +#define FDIMPAR_H #include "cfortran.h" #include "Rtypes.h" @@ -21,6 +21,7 @@ extern "C" { //* mostck = stack dimension for optical photons * //* mxprsn = secondary stack dimension for resonance generator * //* mxpdpm = secondary stack dimension for dpm generators * +//* Mxpscs = secondary stack dimension overall * //* mxoutu = maximum number of output units * //* nallwp = number of allowed particles * //* nelemx = number of maximum allowed elements of a compound * @@ -48,13 +49,14 @@ extern "C" { //*----------------------------------------------------------------------* //* * const Int_t mxxrgn = 10000; -const Int_t mxxmdf = 210; -const Int_t mxxmde = 202; +const Int_t mxxmdf = 510; +const Int_t mxxmde = 502; const Int_t mfstck = 2500; const Int_t mestck = 100; const Int_t mostck = 2000; const Int_t mxprsn = 100; const Int_t mxpdpm = 800; +const Int_t mxpscs = 1999; const Int_t mxoutu = 50; const Int_t nallwp = 64; const Int_t nelemx = 80; diff --git a/TFluka/Femfstk.h b/TFluka/Femfstk.h index cbc6e1173e9..275d59663cf 100644 --- a/TFluka/Femfstk.h +++ b/TFluka/Femfstk.h @@ -38,6 +38,7 @@ typedef struct { Double_t wsnrml[mestck]; Double_t wtemf[mestck]; // weight Double_t agemf[mestck]; // age + Double_t cmpemf[mestck]; Double_t espark[mestck][mkbmx1]; Int_t iespak[mestck][mkbmx2]; Int_t ichemf[mestck]; // charge diff --git a/TFluka/Ffheavy.h b/TFluka/Ffheavy.h index b51e113dfec..5f219e8f684 100644 --- a/TFluka/Ffheavy.h +++ b/TFluka/Ffheavy.h @@ -57,7 +57,7 @@ extern "C" { //*----------------------------------------------------------------------* //* const Int_t mxheav = 100; -const Int_t kxheav = 12; +const Int_t kxheav = 30; typedef struct { Double_t cxheav[mxheav]; diff --git a/TFluka/Ffinuc.h b/TFluka/Ffinuc.h index 6f9501595b5..80fe77b5156 100644 --- a/TFluka/Ffinuc.h +++ b/TFluka/Ffinuc.h @@ -1,3 +1,10 @@ +#ifndef FFINUC +#define FFINUC_H 1 + +#include "Rtypes.h" +#include "cfortran.h" + +#include "Fdimpar.h" extern "C" { //*$ create finuc.add //*copy finuc @@ -52,7 +59,7 @@ extern "C" { //* * //*----------------------------------------------------------------------* //* -const Int_t mxp = 999; +const Int_t mxp = mxpscs; //* typedef struct { @@ -78,3 +85,4 @@ typedef struct { #define FINUC COMMON_BLOCK(FINUC,finuc) COMMON_BLOCK_DEF(finucCommon,FINUC); } +#endif diff --git a/TFluka/Fpaprop.h b/TFluka/Fpaprop.h index 27c5db6ef1e..72940e35109 100644 --- a/TFluka/Fpaprop.h +++ b/TFluka/Fpaprop.h @@ -55,6 +55,7 @@ extern "C" { //*----------------------------------------------------------------------* //* +const Int_t mxgnpr = 33; typedef struct { Double_t am[nallwp+7]; //(-6:NALLWP) Double_t amdisc[nallwp+7]; //(-6:NALLWP) @@ -80,7 +81,7 @@ COMMON_BLOCK_DEF(papropCommon,PAPROP); typedef struct { Char_t btype[nallwp+7][8]; //(-6:NALLWP) - Char_t genpar[30][8]; //(30) + Char_t genpar[mxgnpr][8]; //(30) } chpprpCommon; #define CHPPRP COMMON_BLOCK(CHPPRP,chpprp) COMMON_BLOCK_DEF(chpprpCommon,CHPPRP); diff --git a/TFluka/Fstack.h b/TFluka/Fstack.h index e601507f72d..b3d384d2370 100644 --- a/TFluka/Fstack.h +++ b/TFluka/Fstack.h @@ -72,6 +72,8 @@ extern "C" { //* raddly = delay (s) in production wrt the nominal primary "0" * //* time for particle produced in radioactive decays * //* (i.e. those coming from decays of daughter isotopes) * +//* cmpath = cumulative path travelled by the particle since it * +//* was produced (cm) //* sparek = spare real variables available for k.w.burn * //* ispark = spare integer variables available for k.w.burn * //* ilo = type of the particle (see btype in /paprop/) * @@ -116,6 +118,7 @@ typedef struct { Double_t agestk[mfstck+1]; //(0:MFSTCK) Double_t aknshr[mfstck+1]; //(0:MFSTCK) Double_t raddly[mfstck+1]; //(0:MFSTCK) + Double_t cmpath[mfstck+1]; //(0:MFSTCK) Double_t sparek[mfstck+1][mkbmx1]; //(MKBMX1,0:MFSTCK) Int_t ispark[mfstck+1][mkbmx2]; //(MKBMX2,0:MFSTCK) Int_t ilo[mfstck+1]; //(0:MFSTCK) diff --git a/TFluka/Ftrackr.h b/TFluka/Ftrackr.h index 24b575fe88e..9d7e8451be9 100644 --- a/TFluka/Ftrackr.h +++ b/TFluka/Ftrackr.h @@ -34,6 +34,7 @@ extern "C" { //* ttrack = length of the ith track segment * //* 1 < j < mtrack * //* dtrack = energy deposition of the jth deposition event * +//* dptrck = momentum loss of the jth deposition event * //* * //* jtrack = identity number of the particle * //* etrack = total energy of the particle * @@ -46,6 +47,7 @@ extern "C" { //* biasing techniques are used (for example inelastic * //* interaction length biasing) * //* ctrack = total curved path * +//* cmtrck = cumulative curved path since particle birth * //* zfftrk = of the particle * //* zfrttk = actual z_eff of the particle * //* atrack = age of the particle * @@ -90,6 +92,7 @@ typedef struct { Double_t ztrack[mxtrck+1]; Double_t ttrack[mxtrck]; Double_t dtrack[mxtrck]; + Double_t dptrck[mxtrck][3]; Double_t etrack; Double_t ptrack; Double_t cxtrck; @@ -103,6 +106,7 @@ typedef struct { Double_t zfrttk; Double_t atrack; Double_t ctrack; + Double_t cmtrck; Double_t akshrt; Double_t aklong; Double_t wscrng; diff --git a/TFluka/TFluka.cxx b/TFluka/TFluka.cxx index 18d648d4151..41a32bff4c3 100644 --- a/TFluka/TFluka.cxx +++ b/TFluka/TFluka.cxx @@ -147,7 +147,7 @@ void TFluka::Init() { cout << "\t* InitPhysics() - Prepare input file to be called" << endl; InitPhysics(); // prepare input file - cout << "\t* InitPhysics() - Prepare input file called" << endl; + cout << "\t* InitPhysics() - Prepare input file was called" << endl; if (fVerbosityLevel >=2) cout << "\t* Changing lfdrtr = (" << (GLOBAL.lfdrtr?'T':'F') @@ -470,11 +470,10 @@ Int_t TFluka::IdFromPDG(Int_t pdg) const Int_t TFluka::PDGFromId(Int_t id) const { - // // Return PDG code and pseudo ENDF code from Fluka code - //IPTOKP array goes from official to internal + // IPTOKP array goes from official to internal if (id == -1) { // Cerenkov photon @@ -482,13 +481,13 @@ Int_t TFluka::PDGFromId(Int_t id) const printf("\n PDGFromId: Cerenkov Photon \n"); return 50000050; } - +// Error id if (id == 0) { if (fVerbosityLevel >= 1) printf("PDGFromId: Error id = 0\n"); return -1; } - +// Good id Int_t intfluka = GetFlukaIPTOKP(id); if (intfluka == 0) { if (fVerbosityLevel >= 1) @@ -501,6 +500,7 @@ Int_t TFluka::PDGFromId(Int_t id) const } if (fVerbosityLevel >= 3) printf("mpdgha called with %d %d \n", id, intfluka); + // MPDGHA() goes from fluka internal to pdg. return mpdgha(intfluka); } @@ -1483,6 +1483,7 @@ NOBREM: } } // else if (strncmp(&sProcessFlag[i][0],"PHOT",4) == 0) + // Rayleigh scattering // G3 default value: 0 // G4 process: G4OpRayleigh @@ -1526,15 +1527,83 @@ NOBREM: AliceInp << endl; } } // end of else if (strncmp(&sProcessFlag[i][0],"RAYL",4) == 0) - - else { // processes not yet treated + + // synchrotron radiation in magnetic field + // G3 default value: 0 + // G4 process: G4SynchrotronRadiation + // + // Particles: ?? + // Physics: Not set + // flag = 0 no synchrotron radiation + // flag = 1 synchrotron radiation + //xx gMC ->SetProcess("SYNC",1); // synchrotron radiation generation + else if (strncmp(&sProcessFlag[i][0],"SYNC",4) == 0) { + AliceInp << "*"; + AliceInp << endl; + AliceInp << "*Synchrotron radiation generation is NOT implemented in FLUKA"; + AliceInp << endl; + AliceInp << "*No FLUKA card generated"; + AliceInp << endl; + } + // Automatic calculation of tracking medium parameters // flag = 0 no automatic calculation // flag = 1 automatic calculation //xx gMC ->SetProcess("AUTO",1); // ??? automatic computation of the tracking medium parameters - + else if (strncmp(&sProcessFlag[i][0],"AUTO",4) == 0) { + AliceInp << "*"; + AliceInp << endl; + AliceInp << "*Automatic calculation of tracking medium parameters is always ON in FLUKA"; + AliceInp << endl; + AliceInp << "*No FLUKA card generated"; + AliceInp << endl; + } + + + // To control energy loss fluctuation model + // flag = 0 Urban model + // flag = 1 PAI model + // flag = 2 PAI+ASHO model (not active at the moment) + //xx gMC ->SetProcess("STRA",1); // ??? energy fluctuation model + else if (strncmp(&sProcessFlag[i][0],"STRA",4) == 0) { + if (iProcessValue[i] == 0 || iProcessValue[i] == 2 || iProcessValue[i] == 3) { + AliceInp << "*"; + AliceInp << endl; + AliceInp << "*Ionization energy losses calculation is activated"; + AliceInp << endl; + AliceInp << "*Generated from call: SetProcess('STRA',n);, n=0,1,2"; + AliceInp << endl; + AliceInp << setw(10) << "IONFLUCT "; + AliceInp << setiosflags(ios::fixed) << setiosflags(ios::showpoint) << setprecision(1); + AliceInp << setw(10) << 1.0; // restricted energy loss fluctuations + // (for hadrons and muons) switched on + AliceInp << setw(10) << 1.0; // restricted energy loss fluctuations + // (for e+ and e-) switched on + AliceInp << setw(10) << 1.0; // minimal accuracy + AliceInp << setw(10) << 3.0; // upper bound of the material indices in + // which the respective thresholds apply + AliceInp << setprecision(2); + AliceInp << setw(10) << fLastMaterial; // upper bound of the material indices in which the respective thresholds apply + AliceInp << setprecision(1); + AliceInp << setw(10) << 1.0; // step length in assigning indices + AliceInp << endl; + } + else { + AliceInp << "*"; + AliceInp << endl; + AliceInp << "*Illegal flag value in SetProcess('STRA',?) call."; + AliceInp << endl; + AliceInp << "*No FLUKA card generated"; + AliceInp << endl; + } + } // else if (strncmp(&sProcessFlag[i][0],"STRA",4) == 0) + + + + + else { // processes not yet treated // light photon absorption (Cerenkov photons) // it is turned on when Cerenkov process is turned on @@ -1548,21 +1617,6 @@ NOBREM: // gMC ->SetProcess("LABS",2); // ??? Cerenkov light absorption - // To control energy loss fluctuation model - // flag = 0 Urban model - // flag = 1 PAI model - // flag = 2 PAI+ASHO model (not active at the moment) - //xx gMC ->SetProcess("STRA",1); // ??? energy fluctuation model - - // synchrotron radiation in magnetic field - // G3 default value: 0 - // G4 process: G4SynchrotronRadiation - // - // Particles: ?? - // Physics: Not set - // flag = 0 no synchrotron radiation - // flag = 1 synchrotron radiation - //xx gMC ->SetProcess("SYNC",1); // ??? synchrotron radiation generation cout << "SetProcess for flag=" << &sProcessFlag[i][0] << " value=" << iProcessValue[i] << " not yet implemented!" << endl; } @@ -2039,20 +2093,10 @@ Double_t TFluka::TrackStep() const Double_t TFluka::TrackLength() const { -// Still wrong !!! -// This is the sum of substeps !!! -// TRACKR.ctrack = total curved path of the current step -// Sum of the substeps is identical to TRACKR.ctrack if the is no mag. field -// The sum of all step length starting from the beginning of the track -// for the time being returns only the length in centimeters of the current step - Double_t sum = 0; +// TRACKR.cmtrck = cumulative curved path since particle birth Int_t caller = GetCaller(); - if (caller == 1 || caller == 3 || caller == 4 || caller == 6) { //bxdraw,endraw,mgdraw,usdraw - for ( Int_t j=0;j 0.0) + return 1; + else + return 0; + } else - return 0; // ?????????????? + return 0; } Bool_t TFluka::IsTrackInside() const diff --git a/TFluka/bxdraw.cxx b/TFluka/bxdraw.cxx index c5cdc7feb08..7db50079d87 100644 --- a/TFluka/bxdraw.cxx +++ b/TFluka/bxdraw.cxx @@ -28,7 +28,7 @@ void bxdraw(Int_t& icode, Int_t& mreg, Int_t& newreg, (TVirtualMCApplication::Instance())->Stepping(); fluka->SetCaller(11); fluka->SetTrackIsEntering(); - printf("bxdraw (en) \n"); + printf("bxdraw (en) mreg=%d newreg=%d \n",mreg,newreg); fluka->SetMreg(newreg); (TVirtualMCApplication::Instance())->Stepping(); // fluka->SetCaller(1); diff --git a/TFluka/source_bg.cxx b/TFluka/source_bg.cxx index b80fd600aaa..030aedc9a2f 100644 --- a/TFluka/source_bg.cxx +++ b/TFluka/source_bg.cxx @@ -97,6 +97,8 @@ extern "C" { // projectile gener->SetProjectile("P", 1, 1); gener->SetTarget ("A", 16, 8); + gener->SetBoostLHC(0); + // tell hijing to keep the full parent child chain gener->KeepFullEvent(); // enable jet quenching @@ -119,10 +121,10 @@ extern "C" { gener->Generate(); Int_t npart = stack->GetNprimary(); // Vertex - Float_t za = 4000. * gRandom->Rndm() - 2000.; + Float_t za = 4000. * gRandom->Rndm() -2000.; // Direction - Float_t dir = (gRandom->Rndm() < 0.5) ? 1. : -1.; - +// Float_t dir = (za < 0.) ? 1. : -1.; + Float_t dir = (gRandom->Rndm() < 0.5)? 1. : -1; for (Int_t part=0; partParticle(part); -- 2.31.1