Updates needed for geant4.6 and fluka2004.
authormorsch <morsch@f7af4fe6-9843-0410-8265-dc069ae4e863>
Mon, 26 Jan 2004 16:29:58 +0000 (16:29 +0000)
committermorsch <morsch@f7af4fe6-9843-0410-8265-dc069ae4e863>
Mon, 26 Jan 2004 16:29:58 +0000 (16:29 +0000)
TFluka/Fbeam.h
TFluka/Fdimpar.h
TFluka/Femfstk.h
TFluka/Ffheavy.h
TFluka/Ffinuc.h
TFluka/Fpaprop.h
TFluka/Fstack.h
TFluka/Ftrackr.h
TFluka/TFluka.cxx
TFluka/bxdraw.cxx
TFluka/source_bg.cxx

index 5617bb7..fe94ae1 100644 (file)
@@ -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);
index f728380..98dffaf 100644 (file)
@@ -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;
index cbc6e11..275d596 100644 (file)
@@ -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
index b51e113..5f219e8 100644 (file)
@@ -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];
index 6f95015..80fe77b 100644 (file)
@@ -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
index 27c5db6..72940e3 100644 (file)
@@ -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);
index e601507..b3d384d 100644 (file)
@@ -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)
index 24b575f..9d7e845 100644 (file)
@@ -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 = <z_eff> 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;
index 18d648d..41a32bf 100644 (file)
@@ -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<TRACKR.ntrack;j++) {
-      sum +=TRACKR.ttrack[j];
-    }
-    return sum;
-  }
+  if (caller == 1 || caller == 3 || caller == 4 || caller == 6) //bxdraw,endraw,mgdraw,usdraw
+    return TRACKR.cmtrck;
   else 
     return -1.0;
 }
@@ -2139,14 +2183,16 @@ Double_t TFluka::Etot() const
 //
 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
+// True if the track has positive cummulative length
   Int_t caller = GetCaller();
-  if (caller == 1)
-    return 1; // how to handle double step ?????????????
+  if (caller != 2) { // not eedraw
+    if (TRACKR.cmtrck > 0.0)
+      return 1; 
+    else
+      return 0; 
+  }
   else
-    return 0; // ??????????????
+    return 0;
 }
 
 Bool_t   TFluka::IsTrackInside() const
index c5cdc7f..7db5007 100644 (file)
@@ -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);
index b80fd60..030aedc 100644 (file)
@@ -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; part<npart; part++) {
              particle = stack->Particle(part);