Updates needed for FLUKA2005.6
authormorsch <morsch@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 1 Sep 2005 14:43:05 +0000 (14:43 +0000)
committermorsch <morsch@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 1 Sep 2005 14:43:05 +0000 (14:43 +0000)
TFluka/TFluka.cxx
TFluka/TFlukaConfigOption.cxx
TFluka/TFlukaMCGeometry.cxx
TFluka/abscff.cxx
TFluka/mgdraw.cxx
TFluka/source.cxx
TFluka/stupre.cxx
TFluka/stuprf.cxx

index 2f828ad..536701b 100644 (file)
@@ -34,8 +34,8 @@
 #include "TFluka.h"
 #include "TCallf77.h"      //For the fortran calls
 #include "Fdblprc.h"       //(DBLPRC) fluka common
-#include "Fepisor.h"       //(EPISOR) fluka common
-#include "Ffinuc.h"        //(FINUC)  fluka common
+#include "Fsourcm.h"       //(SOURCM) fluka common
+#include "Fgenstk.h"       //(GENSTK)  fluka common
 #include "Fiounit.h"       //(IOUNIT) fluka common
 #include "Fpaprop.h"       //(PAPROP) fluka common
 #include "Fpart.h"         //(PART)   fluka common
@@ -43,7 +43,7 @@
 #include "Fpaprop.h"       //(PAPROP) fluka common
 #include "Ffheavy.h"       //(FHEAVY) fluka common
 #include "Fopphst.h"       //(OPPHST) fluka common
-#include "Fstack.h"        //(STACK)  fluka common
+#include "Fflkstk.h"       //(FLKSTK) fluka common
 #include "Fstepsz.h"       //(STEPSZ) fluka common
 #include "Fopphst.h"       //(OPPHST) fluka common
 
@@ -307,7 +307,7 @@ void TFluka::ProcessEvent() {
     if (fVerbosityLevel >=3)
        cout << "==> TFluka::ProcessEvent() called." << endl;
     fApplication->GeneratePrimaries();
-    EPISOR.lsouit = true;
+    SOURCM.lsouit = true;
     flukam(1);
     if (fVerbosityLevel >=3)
        cout << "<== TFluka::ProcessEvent() called." << endl;
@@ -598,16 +598,6 @@ void TFluka::Matrix(Int_t& krot, Double_t thetaX, Double_t phiX,
 void TFluka::Gstpar(Int_t itmed, const char* param, Double_t parval) {
 //
 //
-// Check if material is used    
-   if (fVerbosityLevel >= 3) 
-       printf("Gstpar called with %6d %5s %12.4e %6d\n", itmed, param, parval, fGeom->GetFlukaMaterial(itmed));
-   Int_t* reglist;
-   Int_t nreg;
-   reglist = fGeom->GetMaterialList(fGeom->GetFlukaMaterial(itmed), nreg);
-   if (nreg == 0) {
-       return;
-   }
-   
 //
    Bool_t process = kFALSE;
    if (strncmp(param, "DCAY",  4) == 0 ||
@@ -627,10 +617,11 @@ void TFluka::Gstpar(Int_t itmed, const char* param, Double_t parval) {
    {
        process = kTRUE;
    } 
+   
    if (process) {
-       SetProcess(param, Int_t (parval), fGeom->GetFlukaMaterial(itmed));
+       SetProcess(param, Int_t (parval), itmed);
    } else {
-       SetCut(param, parval, fGeom->GetFlukaMaterial(itmed));
+       SetCut(param, parval, itmed);
    }
 }    
 
@@ -1550,11 +1541,11 @@ Int_t TFluka::NSecondaries() const
 
 {
 // Number of secondary particles generated in the current step
-// FINUC.np = number of secondaries except light and heavy ions
+// GENSTK.np = number of secondaries except light and heavy ions
 // FHEAVY.npheav = number of secondaries for light and heavy secondary ions
     Int_t caller = GetCaller();
     if (caller == 6)  // valid only after usdraw
-       return FINUC.np + FHEAVY.npheav;
+       return GENSTK.np + FHEAVY.npheav;
     else if (caller == 50) {
        // Cerenkov Photon production
        return fNCerenkov;
@@ -1571,21 +1562,21 @@ void TFluka::GetSecondary(Int_t isec, Int_t& particleId,
 
     Int_t caller = GetCaller();
     if (caller == 6) {  // valid only after usdraw
-       if (FINUC.np > 0) {
+       if (GENSTK.np > 0) {
            // Hadronic interaction
-           if (isec >= 0 && isec < FINUC.np) {
-               particleId = PDGFromId(FINUC.kpart[isec]);
+           if (isec >= 0 && isec < GENSTK.np) {
+               particleId = PDGFromId(GENSTK.kpart[isec]);
                position.SetX(fXsco);
                position.SetY(fYsco);
                position.SetZ(fZsco);
                position.SetT(TRACKR.atrack);
-               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]);
+               momentum.SetPx(GENSTK.plr[isec]*GENSTK.cxr[isec]);
+               momentum.SetPy(GENSTK.plr[isec]*GENSTK.cyr[isec]);
+               momentum.SetPz(GENSTK.plr[isec]*GENSTK.czr[isec]);
+               momentum.SetE(GENSTK.tki[isec] + PAPROP.am[GENSTK.kpart[isec]+6]);
            }
-           else if (isec >= FINUC.np && isec < FINUC.np + FHEAVY.npheav) {
-               Int_t jsec = isec - FINUC.np;
+           else if (isec >= GENSTK.np && isec < GENSTK.np + FHEAVY.npheav) {
+               Int_t jsec = isec - GENSTK.np;
                particleId = FHEAVY.kheavy[jsec]; // this is Fluka id !!!
                position.SetX(fXsco);
                position.SetY(fYsco);
@@ -1912,7 +1903,7 @@ Double_t TFluka::ParticleLifeTime(Int_t pdg) const
 {
     // Return particle lifetime for particle with pdg code pdg.
     Int_t ifluka = IdFromPDG(pdg);
-    return (PAPROP.thalf[ifluka+6]);
+    return (PAPROP.tmnlf[ifluka+6]);
 }
 
 void TFluka::Gfpart(Int_t pdg, char* name, Int_t& type, Float_t& mass, Float_t& charge, Float_t& tlife)
@@ -1942,30 +1933,30 @@ void TFluka::PrintHeader()
 
 
 
-#define pushcerenkovphoton pushcerenkovphoton_
-#define usersteppingckv    usersteppingckv_
+#define pshckp pshckp_
+#define ustckv ustckv_
 
 
 extern "C" {
-    void pushcerenkovphoton(Double_t & px, Double_t & py, Double_t & pz, Double_t & e,
-                           Double_t & vx, Double_t & vy, Double_t & vz, Double_t & tof,
-                           Double_t & polx, Double_t & poly, Double_t & polz, Double_t & wgt, Int_t& ntr)
-    {
-       //
-       // Pushes one cerenkov photon to the stack
-       //
-       
-       TFluka* fluka =  (TFluka*) gMC;
-       TVirtualMCStack* cppstack = fluka->GetStack();
-       Int_t parent =  TRACKR.ispusr[mkbmx2-1];
-       cppstack->PushTrack(0, parent, 50000050,
-                           px, py, pz, e,
-                            vx, vy, vz, tof,
-                           polx, poly, polz,
-                           kPCerenkov, ntr, wgt, 0); 
-    }
-
-    void usersteppingckv(Int_t & nphot, Int_t & mreg, Double_t & x, Double_t & y, Double_t & z)
+  void pshckp(Double_t & px, Double_t & py, Double_t & pz, Double_t & e,
+             Double_t & vx, Double_t & vy, Double_t & vz, Double_t & tof,
+             Double_t & polx, Double_t & poly, Double_t & polz, Double_t & wgt, Int_t& ntr)
+  {
+    //
+    // Pushes one cerenkov photon to the stack
+    //
+    
+    TFluka* fluka =  (TFluka*) gMC;
+    TVirtualMCStack* cppstack = fluka->GetStack();
+    Int_t parent =  TRACKR.ispusr[mkbmx2-1];
+    cppstack->PushTrack(0, parent, 50000050,
+                       px, py, pz, e,
+                       vx, vy, vz, tof,
+                       polx, poly, polz,
+                       kPCerenkov, ntr, wgt, 0);
+  }
+    
+    void ustckv(Int_t & nphot, Int_t & mreg, Double_t & x, Double_t & y, Double_t & z)
     {
        //
        // Calls stepping in order to signal cerenkov production
@@ -1978,8 +1969,7 @@ extern "C" {
        fluka->SetNCerenkov(nphot);
        fluka->SetCaller(50);
        if (fluka->GetVerbosityLevel() >= 3) 
-       printf("userstepping ckv: %10d %10d %13.3f %13.3f %13.2f %s\n", nphot, mreg, x, y, z, fluka->CurrentVolName());
        (TVirtualMCApplication::Instance())->Stepping();
+       
     }
 }
-
index 61c0355..e64d9da 100644 (file)
@@ -77,8 +77,6 @@ void TFlukaConfigOption::SetCut(const char* flagname, Double_t val)
 
 void TFlukaConfigOption::SetProcess(const char* flagname, Int_t flag)
 {
-    printf("SetProcess %s %5d %5d \n", flagname, flag, fMedium);
-    
     // Set a process flag
     const TString process[15] = 
        {"DCAY", "PAIR", "COMP", "PHOT", "PFIS", "DRAY", "ANNI", "BREM", "MUNU", "CKOV", 
@@ -87,8 +85,6 @@ void TFlukaConfigOption::SetProcess(const char* flagname, Int_t flag)
     Int_t i;
     for (i = 0; i < 15; i++) {
        if (process[i].CompareTo(flagname) == 0) {
-           printf("flag %5d\n", i);
-           
            fProcessFlag[i] = flag;
            if (fMedium == -1) fgDProcessFlag[i] = flag;
            break;
@@ -101,8 +97,19 @@ void TFlukaConfigOption::WriteFlukaInputCards()
     // Write the FLUKA input cards for the set of process flags and cuts
     //
     //
-    if (fMedium > -1) {
+    if (fMedium != -1) {
        TFluka* fluka = (TFluka*) gMC;
+       fMedium = fgGeom->GetFlukaMaterial(fMedium);
+       //
+       // Check if material is actually used
+       Int_t* reglist;
+       Int_t nreg;     
+       reglist = fgGeom->GetMaterialList(fMedium, nreg);
+       if (nreg == 0) {
+           printf("Material not used !\n");
+           return;
+       }
+
        TObjArray *matList = fluka->GetFlukaMaterials();
        Int_t nmaterial =  matList->GetEntriesFast();
        TGeoMaterial* material = 0;
@@ -171,13 +178,13 @@ void TFlukaConfigOption::ProcessDCAY()
 void TFlukaConfigOption::ProcessPAIR()
 {
     // Process PAIR option
-    fprintf(fgFile,"*\n* --- PAIR --- Pair production by gammas, muons and hadrons. Flag = %5d, PPCUTM = %13.4g \n", 
-           fProcessFlag[kPAIR], fCutValue[kPPCUTM]);
+    fprintf(fgFile,"*\n* --- PAIR --- Pair production by gammas, muons and hadrons. Flag = %5d, PPCUTM = %13.4g, PPCUTE = %13.4g\n", 
+           fProcessFlag[kPAIR], fCutValue[kCUTELE], fCutValue[kPPCUTM]);
     //
     // gamma -> e+ e-
     //
     if (fProcessFlag[kPAIR] > 0) {
-       fprintf(fgFile,"EMFCUT    %10.1f%10.1f%10.1f%10.1f%10.1f%10.1fPHOT-THR\n",0., 0., 0.,    fCMatMin, fCMatMax, 1.);
+       fprintf(fgFile,"EMFCUT    %10.1f%10.1f%10.4g%10.1f%10.1f%10.1fPHOT-THR\n",0., 0., 0.0, fCMatMin, fCMatMax, 1.);
     } else {
        fprintf(fgFile,"EMFCUT    %10.1f%10.1f%10.4g%10.1f%10.1f%10.1fPHOT-THR\n",0., 0., 1e10,  fCMatMin, fCMatMax, 1.);
     }
@@ -286,7 +293,7 @@ void TFlukaConfigOption::ProcessPHOT()
     //
 
     if (fProcessFlag[kPHOT] > 0) {
-       fprintf(fgFile,"EMFCUT    %10.1f%10.1f%10.1f%10.1f%10.1f%10.1fPHOT-THR\n",0.   , 0., 0.,  fCMatMin, fCMatMax, 1.);
+       fprintf(fgFile,"EMFCUT    %10.4g%10.4g%10.4g%10.1f%10.1f%10.1fPHOT-THR\n",0., 0., 0.,  fCMatMin, fCMatMax, 1.);
     } else {
        fprintf(fgFile,"EMFCUT    %10.1f%10.4g%10.1f%10.1f%10.1f%10.1fPHOT-THR\n",0., 1.e10, 0.,  fCMatMin, fCMatMax, 1.);
     }
index 166d6db..fef8d90 100644 (file)
@@ -31,7 +31,7 @@
 #include "TGeoManager.h" 
 #include "TGeoVolume.h" 
 #include "TObjString.h"
-#include "Fepisor.h"
+#include "Fsourcm.h"
 
 #ifndef WIN32 
 # define idnrwr idnrwr_
@@ -117,7 +117,7 @@ Int_t gNstep = 0;
 
 ClassImp(TFlukaMCGeometry)
 
-TFlukaMCGeometry* TFlukaMCGeometry::fgInstance=0;
+TFlukaMCGeometry* TFlukaMCGeometry::fgInstance= NULL;
 
 //_____________________________________________________________________________
 TFlukaMCGeometry::TFlukaMCGeometry(const char *name, const char *title) 
@@ -128,9 +128,12 @@ TFlukaMCGeometry::TFlukaMCGeometry(const char *name, const char *title)
   //
   fDebug        = kFALSE;
   fLastMaterial = 0;
+  fCurrentRegion   = 0;
+  fCurrentLattice  = 0;
   fNextRegion   = 0;
   fNextLattice  = 0;
   fRegionList   = 0;
+  fIndmat = 0;
   gFluka = (TFluka*)gMC;
   gMCGeom = this;
   gNstep = 0;
@@ -146,9 +149,12 @@ TFlukaMCGeometry::TFlukaMCGeometry()
   //
   fDebug        = kFALSE;
   fLastMaterial = 0;
+  fCurrentRegion   = 0;
+  fCurrentLattice  = 0;
   fNextRegion   = 0;
   fNextLattice  = 0;
   fRegionList   = 0;
+  fIndmat = 0;
   gFluka = (TFluka*)gMC;
   gMCGeom = this;
   gNstep = 0;
@@ -219,11 +225,13 @@ Int_t TFlukaMCGeometry::GetMedium() const
 Int_t TFlukaMCGeometry::GetFlukaMaterial(Int_t imed) const
 {
 // Returns FLUKA material index for medium IMED
-   TGeoMedium *med = (TGeoMedium*)gGeoManager->GetListOfMedia()->At(imed-1);
+     TGeoMedium *med = (TGeoMedium*)gGeoManager->GetListOfMedia()->At(imed-1);
    if (!med) {
       Error("GetFlukaMaterial", "MEDIUM %i nor found", imed);
       return -1;
    }
+   TGeoMaterial* mat = med->GetMaterial();
+   if (!mat->IsUsed()) return -1;
    Int_t imatfl = med->GetMaterial()->GetIndex();
    return imatfl;   
 }
@@ -296,10 +304,10 @@ TGeoMaterial * TFlukaMCGeometry::GetMakeWrongMaterial(Double_t z)
    }
    TGeoMixture *mix = 0;
    TGeoElement *element;
-   TGeoElementTable *table = TGeoElementTable::Instance();
+   TGeoElementTable *table = gGeoManager->GetElementTable();
    switch (ind) {
       case 0: // AIR
-         mix = new TGeoMixture("AIR", 4, 0.001205);
+         mix = new TGeoMixture("TPC_AIR", 4, 0.001205);
          element = table->GetElement(6); // C
          mix->DefineElement(0, element, 0.000124);
          element = table->GetElement(7); // N
@@ -310,7 +318,7 @@ TGeoMaterial * TFlukaMCGeometry::GetMakeWrongMaterial(Double_t z)
          mix->DefineElement(3, element, 0.012827);
          break;
       case 1: //SDD SI CHIP
-         mix = new TGeoMixture("SDD_SI", 6, 2.4485);
+         mix = new TGeoMixture("ITS_SDD_SI", 6, 2.4485);
          element = table->GetElement(1);
          mix->DefineElement(0, element, 0.004367771);
          element = table->GetElement(6);
@@ -325,14 +333,14 @@ TGeoMaterial * TFlukaMCGeometry::GetMakeWrongMaterial(Double_t z)
          mix->DefineElement(5, element, 0.09814344903);
          break;
       case 2:  // WATER
-         mix = new TGeoMixture("WATER", 2, 1.0);
+         mix = new TGeoMixture("ITS_WATER", 2, 1.0);
          element = table->GetElement(1);
          mix->DefineElement(0, element, 0.111898344);
          element = table->GetElement(8);
          mix->DefineElement(1, element, 0.888101656);
          break;
       case 3: // CERAMICS
-         mix = new TGeoMixture("CERAMICS", 5, 3.6);
+         mix = new TGeoMixture("ITS_CERAMICS", 5, 3.6);
          element = table->GetElement(8);
          mix->DefineElement(0, element, 0.59956);
          element = table->GetElement(13);
@@ -345,7 +353,7 @@ TGeoMaterial * TFlukaMCGeometry::GetMakeWrongMaterial(Double_t z)
          mix->DefineElement(4, element, 0.0115);
          break;
       case 4: // EPOXY
-         mix = new TGeoMixture("G10FR4", 4, 1.8);
+         mix = new TGeoMixture("MUON_G10FR4", 4, 1.8);
          element = table->GetElement(1);
          mix->DefineElement(0, element, 0.19);
          element = table->GetElement(6);
@@ -367,7 +375,7 @@ TGeoMaterial * TFlukaMCGeometry::GetMakeWrongMaterial(Double_t z)
          mix->DefineElement(3, element, 0.28);
          break;
       case 6: // KAPTON
-         mix = new TGeoMixture("KAPTON", 4, 1.3);
+         mix = new TGeoMixture("ITS_KAPTON", 4, 1.3);
          element = table->GetElement(1);
          mix->DefineElement(0, element, 0.026363415);
          element = table->GetElement(6);
@@ -378,7 +386,7 @@ TGeoMaterial * TFlukaMCGeometry::GetMakeWrongMaterial(Double_t z)
          mix->DefineElement(3, element, 0.209238060);
          break;
       case 7: // INOX
-         mix = new TGeoMixture("INOX", 9, 7.9);
+         mix = new TGeoMixture("ITS_INOX", 9, 7.9);
          element = table->GetElement(6);
          mix->DefineElement(0, element, 0.0003);
          element = table->GetElement(14);
@@ -410,7 +418,7 @@ TGeoMaterial * TFlukaMCGeometry::GetMakeWrongMaterial(Double_t z)
          mix->DefineElement(3, element, 0.19137459);
          break;
       case 9: // SDD-C-AL
-         mix = new TGeoMixture("SDD-C-AL", 5, 1.9837);
+         mix = new TGeoMixture("ITS_SDD-C-AL", 5, 1.9837);
          element = table->GetElement(1);
          mix->DefineElement(0, element, 0.022632);
          element = table->GetElement(6);
@@ -423,7 +431,7 @@ TGeoMaterial * TFlukaMCGeometry::GetMakeWrongMaterial(Double_t z)
          mix->DefineElement(4, element, 0.1);
          break;
       case 10: // X7R-CAP
-         mix = new TGeoMixture("X7R-CAP", 7, 6.72);
+         mix = new TGeoMixture("ITS_X7R-CAP", 7, 6.72);
          element = table->GetElement(8);
          mix->DefineElement(0, element, 0.085975822);
          element = table->GetElement(22);
@@ -440,14 +448,14 @@ TGeoMaterial * TFlukaMCGeometry::GetMakeWrongMaterial(Double_t z)
          mix->DefineElement(6, element, 0.2081768);
          break;
       case 11: // SDD ruby sph. Al2O3
-         mix = new TGeoMixture("AL2O3", 2, 3.97);
+         mix = new TGeoMixture("ITS_AL2O3", 2, 3.97);
          element = table->GetElement(8);
          mix->DefineElement(0, element, 0.5293);
          element = table->GetElement(13);
          mix->DefineElement(1, element, 0.4707);
          break;
       case 12: // SDD HV microcable
-         mix = new TGeoMixture("HV-CABLE", 5, 1.6087);
+         mix = new TGeoMixture("ITS_HV-CABLE", 5, 1.6087);
          element = table->GetElement(1);
          mix->DefineElement(0, element, 0.01983871336);
          element = table->GetElement(6);
@@ -460,7 +468,7 @@ TGeoMaterial * TFlukaMCGeometry::GetMakeWrongMaterial(Double_t z)
          mix->DefineElement(4, element, 0.247536);
          break;
       case 13: //SDD LV+signal cable
-         mix = new TGeoMixture("LV-CABLE", 5, 2.1035);
+         mix = new TGeoMixture("ITS_LV-CABLE", 5, 2.1035);
          element = table->GetElement(1);
          mix->DefineElement(0, element, 0.0082859922);
          element = table->GetElement(6);
@@ -473,7 +481,7 @@ TGeoMaterial * TFlukaMCGeometry::GetMakeWrongMaterial(Double_t z)
          mix->DefineElement(4, element, 0.68572);
          break;
       case 14: //SDD hybrid microcab
-         mix = new TGeoMixture("HYB-CAB", 5, 2.0502);
+         mix = new TGeoMixture("ITS_HYB-CAB", 5, 2.0502);
          element = table->GetElement(1);
          mix->DefineElement(0, element, 0.00926228815);
          element = table->GetElement(6);
@@ -486,7 +494,7 @@ TGeoMaterial * TFlukaMCGeometry::GetMakeWrongMaterial(Double_t z)
          mix->DefineElement(4, element, 0.64869);
          break;
       case 15: //SDD anode microcab
-         mix = new TGeoMixture("ANOD-CAB", 5, 1.7854);
+         mix = new TGeoMixture("ITS_ANOD-CAB", 5, 1.7854);
          element = table->GetElement(1);
          mix->DefineElement(0, element, 0.0128595919215);
          element = table->GetElement(6);
@@ -499,7 +507,7 @@ TGeoMaterial * TFlukaMCGeometry::GetMakeWrongMaterial(Double_t z)
          mix->DefineElement(4, element, 0.431909);
          break;
       case 16: // inox/alum
-         mix = new TGeoMixture("INOX-AL", 5, 3.0705);
+         mix = new TGeoMixture("ITS_INOX-AL", 5, 3.0705);
          element = table->GetElement(13);
          mix->DefineElement(0, element, 0.816164);
          element = table->GetElement(14);
@@ -511,7 +519,7 @@ TGeoMaterial * TFlukaMCGeometry::GetMakeWrongMaterial(Double_t z)
          element = table->GetElement(28);
          mix->DefineElement(4, element, 0.0183836);
       case 17: // MYLAR
-         mix = new TGeoMixture("MYLAR", 3, 1.39);
+         mix = new TGeoMixture("TPC_MYLAR", 3, 1.39);
          element = table->GetElement(1);
          mix->DefineElement(0, element, 0.0416667);
          element = table->GetElement(6);
@@ -520,31 +528,31 @@ TGeoMaterial * TFlukaMCGeometry::GetMakeWrongMaterial(Double_t z)
          mix->DefineElement(2, element, 0.333333);
          break;
       case 18: // SPDBUS(AL+KPT+EPOX)   - unknown composition
-         mix = new TGeoMixture("SPDBUS", 1, 1.906);
+         mix = new TGeoMixture("ITS_SPDBUS", 1, 1.906);
          element = table->GetElement(9);
          mix->DefineElement(0, element, 1.);
          z = element->Z();
          break;
       case 19: // SDD/SSD rings   - unknown composition
-         mix = new TGeoMixture("SDDRINGS", 1, 1.8097);
+         mix = new TGeoMixture("ITS_SDDRINGS", 1, 1.8097);
          element = table->GetElement(6);
          mix->DefineElement(0, element, 1.);
          z = element->Z();
          break;
       case 20: // SPD end ladder   - unknown composition
-         mix = new TGeoMixture("SPDEL", 1, 3.6374);
+         mix = new TGeoMixture("ITS_SPDEL", 1, 3.6374);
          element = table->GetElement(22);
          mix->DefineElement(0, element, 1.);
          z = element->Z();
          break;
       case 21: // SDD end ladder   - unknown composition
-         mix = new TGeoMixture("SDDEL", 1, 0.3824);
+         mix = new TGeoMixture("ITS_SDDEL", 1, 0.3824);
          element = table->GetElement(30);
          mix->DefineElement(0, element, 1.);
          z = element->Z();
          break;
       case 22: // SSD end ladder   - unknown composition
-         mix = new TGeoMixture("SSDEL", 1, 0.68);
+         mix = new TGeoMixture("ITS_SSDEL", 1, 0.68);
          element = table->GetElement(16);
          mix->DefineElement(0, element, 1.);
          z = element->Z();
@@ -582,7 +590,7 @@ void TFlukaMCGeometry::CreateFlukaMatFile(const char *fname)
    Int_t i,j,idmat;
    Int_t counttothree, nelem;
    Double_t a,z,rho, w;
-   TGeoElementTable *table = TGeoElementTable::Instance();
+   TGeoElementTable *table = gGeoManager->GetElementTable();
    TGeoElement *element;
    element = table->GetElement(13);
    element->SetTitle("ALUMINUM"); // this is how FLUKA likes it ...
@@ -624,7 +632,7 @@ void TFlukaMCGeometry::CreateFlukaMatFile(const char *fname)
    // Adjust material names and add them to FLUKA list
    for (i=0; i<nmater; i++) {
       mat = (TGeoMaterial*)matlist->At(i);
-      if (!mat->IsUsed()) continue;      
+      if (!mat->IsUsed()) continue;
       z = mat->GetZ();
       a = mat->GetA();
       rho = mat->GetDensity();
@@ -634,27 +642,7 @@ void TFlukaMCGeometry::CreateFlukaMatFile(const char *fname)
       } 
       matname = mat->GetName();
       FlukaMatName(matname);
-/*
-      // material with one element: create it as mixture since it can be duplicated    
-      if (!mat->IsMixture()) {
-         // normal materials
-         mix = new TGeoMixture(matname.Data(), 1, rho);
-         mix->DefineElement(0, mat->GetElement(), 1.);
-         mat->SetIndex(nfmater+3);
-         for (j=0; j<nmed; j++) {
-            med = (TGeoMedium*)medlist->At(j);
-            if (med->GetMaterial() == mat) {
-               med->SetMaterial(mix);
-               if (mat->GetCerenkovProperties()) {
-                  mix->SetCerenkovProperties(mat->GetCerenkovProperties());
-                  mat->SetCerenkovProperties(0);
-               }
-               break;
-            }
-         }                              
-         mat = (TGeoMaterial*)mix;
-      }
-*/    
+
       mat->SetIndex(nfmater+3);
       objstr = new TObjString(matname.Data());
       fMatList->Add(mat);
@@ -800,11 +788,6 @@ void TFlukaMCGeometry::CreateFlukaMatFile(const char *fname)
    }
    out.close();
    fLastMaterial = nfmater+2;
-   
-   if (!gFluka->IsGeneratePemf()) {
-       if (gSystem->AccessPathName("FlukaVmc.pemf")) Fatal("CreateFlukaMatFile", "No pemf file in working directory");
-       return;
-   }
 }
 
 void TFlukaMCGeometry::CreatePemfFile()
@@ -1278,6 +1261,18 @@ void TFlukaMCGeometry::ToFlukaString(TString &str) const
 //_____________________________________________________________________________
 void TFlukaMCGeometry::FlukaMatName(TString &str) const
 {
+// Strip the detector name
+     
+    TObjArray * tokens = str.Tokenize("_");
+    Int_t ntok = tokens->GetEntries();
+    if (ntok > 0) {
+       TString head = ((TObjString*) tokens->At(0))->GetString();
+       Int_t nhead = head.Length();
+       str = str.Remove(0, nhead + 1);
+    }
+    tokens->Clear();
+    delete tokens;
+
 // Convert a name to upper case 8 chars.
    ToFlukaString(str);
    Int_t ilast;
@@ -1322,6 +1317,7 @@ void TFlukaMCGeometry::Vname(const char *name, char *vname) const
 }
 
 
+
 // FLUKA GEOMETRY WRAPPERS - to replace FLUGG wrappers
 
 //_____________________________________________________________________________
@@ -1361,10 +1357,16 @@ void g1wr(Double_t &pSx, Double_t &pSy, Double_t &pSz,
       printf("   oldReg=%i  oldLttc=%i  pstep=%f\n",oldReg, oldLttc, propStep);
    }
    Int_t olttc = oldLttc;
-   if (oldLttc<0) {
+   if (oldLttc<=0) {
       gGeoManager->FindNode(pSx,pSy,pSz);
       olttc = gGeoManager->GetCurrentNodeId()+1;
+      if (gMCGeom->IsDebugging()) {
+         printf("WOOPS: old reg/latt = %i/%i\n",oldReg,oldLttc);
+         printf("point: (%16.12f, %16.12f, %16.12f) in lttc=%i\n", pSx,pSy,pSz,olttc);
+      }   
    }   
+   Int_t ccreg,cclat;
+   gMCGeom->GetCurrentRegion(ccreg,cclat);
    gMCGeom->SetCurrentRegion(oldReg, olttc);
    // Initialize default return values
    lttcFlag = 0;
@@ -1417,31 +1419,46 @@ void g1wr(Double_t &pSx, Double_t &pSy, Double_t &pSz,
    Double_t extra = 1.E-10;
 
      
-//      printf("ERROR: (%f, %f, %f)\n",pSx,pSy,pSz);
    if (gMCGeom->IsDebugging()) printf("   current path: %s\n", gGeoManager->GetPath());
    gGeoManager->SetCurrentPoint(pSx+extra*pV[0], pSy+extra*pV[1], pSz+extra*pV[2]);
    gGeoManager->SetCurrentDirection(pV);
-   gGeoManager->FindNextBoundary(-propStep);
+//   gGeoManager->FindNextBoundary(-propStep);
+   gGeoManager->FindNextBoundary(-100000.);
    Double_t snext = gGeoManager->GetStep();
 
    if (snext<=0.0) {
-      saf = 0.0;
-      newReg = -3;
-      sLt[lttcFlag] = 0.0;
-      if (gMCGeom->IsDebugging()) printf("BACK SCATTERING\n");
-      return;
+      snext = 0.0;
+      /*
+      if (!crossed) {
+         // artefact due to MS
+         saf = 0.0;
+         newReg = -3;
+         sLt[lttcFlag] = 0.0;
+         if (gMCGeom->IsDebugging()) {
+            printf("BACK SCATTERING (reg/lttc = %i/%i)\n", oldReg,oldLttc);
+            printf("point: (%16.11f, %16.11f, %16.11f)\n", pSx,pSy,pSz);
+            printf("dir :  (%16.11f, %16.11f, %16.11f)\n", pV[0],pV[1],pV[2]);
+         }   
+         return;
+      }
+      */
+      // else ... 
+      // FLUKA native detects some extra errors even when a boundary was crossed, returning a -3
+      // If the particle turns back to the original region before crossing at the first step,
+      // we just return the distance to the boundary, not issuing an error (not due to MS)
    }  
-
    snext += extra;
    saf = gGeoManager->GetSafeDistance();
    saf -= extra;
+   if (saf>snext) printf("ERROR: saf=%f .GT. snext=%f\n", saf,snext); 
    if (saf<0) saf=0.0;
    else       saf -= saf*3.0e-09;
+//   saf *= 0.3;
+   PAREM.dist = snext;
    NORLAT.distn = snext;
    NORLAT.xn[0] += snext*pV[0];
    NORLAT.xn[1] += snext*pV[1];
    NORLAT.xn[2] += snext*pV[2];
-//   saf = 0.0; // !!! TEMPORARY FOR TESTING MAGSPHF - TO BE REMOVED 
    if (snext>propStep) {
    // Next boundary further than proposed step, which is approved
       retStep = propStep;
@@ -1449,6 +1466,8 @@ void g1wr(Double_t &pSx, Double_t &pSy, Double_t &pSz,
       return;
    }
    // The next boundary is closer. We try to cross it.
+   if (saf>propStep) saf = propStep; // Safety should be less than the proposed step if a boundary will be crossed
+//   saf = 0.0; // !!! SAFETY SHOULD BE 0 IF THE BOUNDARY WILL BE CROSSED ???
    gGeoManager->SetCurrentPoint(pSx,pSy,pSz);
    Double_t *point = gGeoManager->GetCurrentPoint();
    Double_t *dir = gGeoManager->GetCurrentDirection();
@@ -1456,14 +1475,23 @@ void g1wr(Double_t &pSx, Double_t &pSy, Double_t &pSz,
    memcpy(pt, point, 3*sizeof(Double_t));
    
    Int_t i;
-   for (i=0;i<3;i++) point[i] += (snext+1E-6)*dir[i];
+   extra = 1.E-13;
+   for (i=0;i<3;i++) point[i] += (snext+extra)*dir[i];
    // locate next region
    gGeoManager->FindNode();
    newLttc = (gGeoManager->IsOutside())?(TFlukaMCGeometry::kLttcOutside):gGeoManager->GetCurrentNodeId()+1;
+   while (newLttc==olttc) {
+      extra *= 10.;
+      if (extra>1.E-5) break;
+      for (i=0;i<3;i++) point[i] += extra*dir[i];
+      gGeoManager->FindNode();
+      newLttc = (gGeoManager->IsOutside())?(TFlukaMCGeometry::kLttcOutside):gGeoManager->GetCurrentNodeId()+1;
+   }           
    gGeoManager->SetCurrentPoint(pt);
    newReg = (gGeoManager->IsOutside())?(gMCGeom->NofVolumes()+1):gGeoManager->GetCurrentVolume()->GetNumber();
    if (gMCGeom->IsDebugging()) printf("   newReg=%i newLttc=%i\n", newReg, newLttc);
 
+
    // We really crossed the boundary, but is it the same region ?
    gMCGeom->SetNextRegion(newReg, newLttc);
    if (newReg==oldReg && newLttc!=olttc) {
@@ -1583,12 +1611,14 @@ void lkwr(Double_t &pSx, Double_t &pSy, Double_t &pSz,
       printf("   in: pV=(%f, %f, %f)\n", pV[0], pV[1], pV[2]);
       printf("   in: oldReg=%i oldLttc=%i\n", oldReg, oldLttc);
    }   
+/*
    NORLAT.xn[0] = pSx;
    NORLAT.xn[1] = pSy;
    NORLAT.xn[2] = pSz;
    NORLAT.wn[0] = pV[0];
    NORLAT.wn[1] = pV[1];
    NORLAT.wn[2] = pV[2];
+*/
    TGeoNode *node = gGeoManager->FindNode(pSx, pSy, pSz);
    if (gGeoManager->IsOutside()) {
       newReg = gMCGeom->NofVolumes()+1;
index eabb2e5..94967ba 100644 (file)
@@ -14,6 +14,8 @@
 extern "C" {
 Double_t abscff(Double_t& wvlngt, Double_t& /*omgpho*/, Int_t& mmat)
 {
+//    printf("abscff%f %d\n", wvlngt, mmat);
+    
 //
 //  Return absorption length  for given photon energy and material
 //
index 3dbe4f3..f2ef1cb 100644 (file)
@@ -9,7 +9,7 @@
 #include "Fdblprc.h"  //(DBLPRC) fluka common
 #include "Ftrackr.h"  //(TRACKR) fluka common
 #include "Fopphst.h"  //(OPPHST) fluka common
-#include "Fstack.h"   //(STACK)  fluka common
+#include "Fflkstk.h"  //(FLKSTK) fluka common
 
 #ifndef WIN32
 # define mgdraw mgdraw_
@@ -31,7 +31,7 @@ void mgdraw(Int_t& icode, Int_t& mreg)
     if (TRACKR.jtrack == -1) {
        trackId = OPPHST.louopp[OPPHST.lstopp];
        if (trackId == 0) {
-           trackId = STACK.ispark[STACK.lstack][mkbmx2-1];
+           trackId = FLKSTK.ispark[FLKSTK.npflka][mkbmx2-1];
        }
     } else {
        trackId = TRACKR.ispusr[mkbmx2-1];
@@ -48,10 +48,14 @@ void mgdraw(Int_t& icode, Int_t& mreg)
     if (!TRACKR.ispusr[mkbmx2 - 2]) {
        //
        // Single step
-       if (verbosityLevel >= 3) {
+       if (TRACKR.jtrack == -1 && trackId == 109340) {
            cout << endl << " !!! I am in mgdraw - calling Stepping(): " << icode << endl;
            cout << endl << " Track Id = " << trackId << " region = " << mreg << endl;
+           printf("Stepsize %13.5e \n", fluka->TrackStep());
        }
+
+
+
       
        (TVirtualMCApplication::Instance())->Stepping();
        fluka->SetTrackIsNew(kFALSE);
index 1a2edcb..f66b562 100644 (file)
@@ -4,14 +4,12 @@
 // Fluka commons
 #include "Fdblprc.h"  //(DBLPRC) fluka common
 #include "Fdimpar.h"  //(DIMPAR) fluka parameters
-#include "Fepisor.h"  //(EPISOR) fluka common
-#include "Fstack.h"   //(STACK)  fluka common
-#include "Fstars.h"   //(STARS)  fluka common
-#include "Fbeam.h"    //(BEAM)   fluka common
+#include "Fsourcm.h"  //(EPISOR) fluka common
+#include "Fflkstk.h"  //(FLKSTK)  fluka common
+#include "Fsumcou.h"   //(SUMCOU)  fluka common
 #include "Fpaprop.h"  //(PAPROP) fluka common
 #include "Fltclcm.h"  //(LTCLCM) fluka common
 #include "Fopphst.h"  //(OPPHST) fluka common
-//#include "Fcaslim.h"  //(CASLIM) fluka common
 
 //Virutal MC
 #include "TFluka.h"
@@ -73,7 +71,7 @@ extern "C" {
     Bool_t debug = (verbosityLevel>=3)?kTRUE:kFALSE;
     if (debug) {
       cout << "==> source(" << nomore << ")" << endl;
-      cout << "\t* EPISOR.lsouit = " << (EPISOR.lsouit?'T':'F') << endl;
+      cout << "\t* SOURCM.lsouit = " << (SOURCM.lsouit?'T':'F') << endl;
     }  
 
     static Bool_t lfirst = true;
@@ -102,9 +100,9 @@ extern "C" {
     }
 
     if (lfirst) {
-       EPISOR.tkesum = zerzer;
+       SOURCM.tkesum = zerzer;
        lfirst = false;
-       EPISOR.lussrc = true;
+       SOURCM.lussrc = true;
     } else {
 //
 // Post-track actions for primary track
@@ -120,9 +118,9 @@ extern "C" {
 
     if (itrack<0) {
       nomore = 1;
-      EPISOR.lsouit = false;
+      SOURCM.lsouit = false;
       if (debug) {
-         cout << "\t* EPISOR.lsouit = " << (EPISOR.lsouit?'T':'F') << endl;
+         cout << "\t* SOURCM.lsouit = " << (SOURCM.lsouit?'T':'F') << endl;
          cout << "\t* No more particles. Exiting..." << endl;
          cout << "<== source(" << nomore << ")" << endl;
       }   
@@ -135,7 +133,7 @@ extern "C" {
        printf("Event has been stopped by user !");
        fluka->SetStopEvent(kFALSE);
        nomore = 1;
-       EPISOR.lsouit = false;
+       SOURCM.lsouit = false;
        return;
     }
 
@@ -158,7 +156,7 @@ extern "C" {
            << particle->Energy() << " GeV "  
            << particle->GetMass() << " GeV " << endl;
     }   
-    /* Lstack is the stack counter: of course any time source is called it
+    /* Npflka is the stack counter: of course any time source is called it
      * must be =0
      */
     /* Cosines (tx,ty,tz)*/
@@ -167,106 +165,105 @@ extern "C" {
     Double_t cosz = TMath::Sqrt(oneone - cosx*cosx - cosy*cosy);
     if (particle->Pz() < 0.) cosz = -cosz;    
 
-    //STACK.ilo[STACK.lstack] = BEAM.ijbeam;
     if (pdg != 50000050 &&  pdg !=  50000051) {
-       STACK.lstack++;
+       FLKSTK.npflka++;
        Int_t ifl =  fluka-> IdFromPDG(pdg);
-       STACK.ilo[STACK.lstack] = ifl;
-       /* Wt is the weight of the particle*/
-       STACK.wt[STACK.lstack] = oneone;
-       STARS.weipri += STACK.wt[STACK.lstack];
+       FLKSTK.iloflk[FLKSTK.npflka] = ifl;
+       /* Wtflk is the weight of the particle*/
+       FLKSTK.wtflk[FLKSTK.npflka] = oneone;
+       SUMCOU.weipri += FLKSTK.wtflk[FLKSTK.npflka];
        
-       STACK.lo[STACK.lstack] = 1;
+       FLKSTK.loflk[FLKSTK.npflka] = 1;
        
        /* User dependent flag:*/
-       STACK.louse[STACK.lstack] = 0;
+       FLKSTK.louse[FLKSTK.npflka] = 0;
 
        /* User dependent spare variables:*/
        Int_t ispr = 0;
        for (ispr = 0; ispr < mkbmx1; ispr++)
-           STACK.sparek[STACK.lstack][ispr] = zerzer;
+           FLKSTK.sparek[FLKSTK.npflka][ispr] = zerzer;
 
        /* User dependent spare flags:*/
        for (ispr = 0; ispr < mkbmx2; ispr++)
-           STACK.ispark[STACK.lstack][ispr] = 0;
+           FLKSTK.ispark[FLKSTK.npflka][ispr] = 0;
 
        /* Save the track number of the stack particle:*/
-       STACK.ispark[STACK.lstack][mkbmx2-1] = itrack;
-       STACK.nparma++;
-       STACK.numpar[STACK.lstack] = STACK.nparma;
-       STACK.nevent[STACK.lstack] = 0;
-       STACK.dfnear[STACK.lstack] = +zerzer;
+       FLKSTK.ispark[FLKSTK.npflka][mkbmx2-1] = itrack;
+       FLKSTK.nparma++;
+       FLKSTK.numpar[FLKSTK.npflka] = FLKSTK.nparma;
+       FLKSTK.nevent[FLKSTK.npflka] = 0;
+       FLKSTK.dfnear[FLKSTK.npflka] = +zerzer;
        
        /* Particle age (s)*/
-       STACK.agestk[STACK.lstack] = +zerzer;
-       STACK.cmpath[STACK.lstack] = +zerzer;
-       STACK.aknshr[STACK.lstack] = -twotwo;
+       FLKSTK.agestk[FLKSTK.npflka] = +zerzer;
+       FLKSTK.cmpath[FLKSTK.npflka] = +zerzer;
+       FLKSTK.aknshr[FLKSTK.npflka] = -twotwo;
 
        /* Group number for "low" energy neutrons, set to 0 anyway*/
-       STACK.igroup[STACK.lstack] = 0;
+       FLKSTK.igroup[FLKSTK.npflka] = 0;
        
        /* Kinetic energy */
        Double_t p    = particle->P();
        Double_t mass = PAPROP.am[ifl + 6];
-       STACK.tke[STACK.lstack] =  TMath::Sqrt( p * p + mass * mass) - mass;
+       FLKSTK.tkeflk[FLKSTK.npflka] =  TMath::Sqrt( p * p + mass * mass) - mass;
        /* Particle momentum*/
-       STACK.pmom [STACK.lstack] = p;
+       FLKSTK.pmoflk [FLKSTK.npflka] = p;
 
-       STACK.tx [STACK.lstack] = cosx;
-       STACK.ty [STACK.lstack] = cosy;
-       STACK.tz [STACK.lstack] = cosz;
+       FLKSTK.txflk [FLKSTK.npflka] = cosx;
+       FLKSTK.tyflk [FLKSTK.npflka] = cosy;
+       FLKSTK.tzflk [FLKSTK.npflka] = cosz;
     
        /* Polarization cosines:*/
        if (polarisation.Mag()) {
            Double_t cospolx = polarisation.Px() / polarisation.Mag();
            Double_t cospoly = polarisation.Py() / polarisation.Mag();
            Double_t cospolz = sqrt(oneone - cospolx * cospolx - cospoly * cospoly);
-           STACK.txpol [STACK.lstack] = cospolx;
-           STACK.typol [STACK.lstack] = cospoly;
-           STACK.tzpol [STACK.lstack] = cospolz;
+           FLKSTK.txpol [FLKSTK.npflka] = cospolx;
+           FLKSTK.typol [FLKSTK.npflka] = cospoly;
+           FLKSTK.tzpol [FLKSTK.npflka] = cospolz;
        }
        else {
-           STACK.txpol [STACK.lstack] = -twotwo;
-           STACK.typol [STACK.lstack] = +zerzer;
-           STACK.tzpol [STACK.lstack] = +zerzer;
+           FLKSTK.txpol [FLKSTK.npflka] = -twotwo;
+           FLKSTK.typol [FLKSTK.npflka] = +zerzer;
+           FLKSTK.tzpol [FLKSTK.npflka] = +zerzer;
        }
        
        /* Particle coordinates*/
        // Vertext coordinates;
-       STACK.xa [STACK.lstack] = particle->Vx();
-       STACK.ya [STACK.lstack] = particle->Vy();
-       STACK.za [STACK.lstack] = particle->Vz();
+       FLKSTK.xflk [FLKSTK.npflka] = particle->Vx();
+       FLKSTK.yflk [FLKSTK.npflka] = particle->Vy();
+       FLKSTK.zflk [FLKSTK.npflka] = particle->Vz();
     
        /*  Calculate the total kinetic energy of the primaries: don't change*/
-       Int_t st_ilo =  STACK.ilo[STACK.lstack];
+       Int_t st_ilo =  FLKSTK.iloflk[FLKSTK.npflka];
        if ( st_ilo != 0 )
-           EPISOR.tkesum += 
-               ((STACK.tke[STACK.lstack] + PAPROP.amdisc[st_ilo+6])
-                * STACK.wt[STACK.lstack]);
+           SOURCM.tkesum += 
+               ((FLKSTK.tkeflk[FLKSTK.npflka] + PAPROP.amdisc[st_ilo+6])
+                * FLKSTK.wtflk[FLKSTK.npflka]);
        else
-           EPISOR.tkesum += (STACK.tke[STACK.lstack] * STACK.wt[STACK.lstack]);
+           SOURCM.tkesum += (FLKSTK.tkeflk[FLKSTK.npflka] * FLKSTK.wtflk[FLKSTK.npflka]);
        
        /*  Here we ask for the region number of the hitting point.
-        *     NREG (LSTACK) = ...
+        *     NRGFLK (LFLKSTK) = ...
         *  The following line makes the starting region search much more
         *  robust if particles are starting very close to a boundary:
         */
-       geocrs( STACK.tx[STACK.lstack], 
-               STACK.ty[STACK.lstack], 
-               STACK.tz[STACK.lstack] );
+       geocrs( FLKSTK.txflk[FLKSTK.npflka], 
+               FLKSTK.tyflk[FLKSTK.npflka], 
+               FLKSTK.tzflk[FLKSTK.npflka] );
     
        Int_t idisc;
 
-       georeg ( STACK.xa[STACK.lstack], 
-                STACK.ya[STACK.lstack], 
-                STACK.za[STACK.lstack],
-                STACK.nreg[STACK.lstack], 
+       georeg ( FLKSTK.xflk[FLKSTK.npflka], 
+                FLKSTK.yflk[FLKSTK.npflka], 
+                FLKSTK.zflk[FLKSTK.npflka],
+                FLKSTK.nrgflk[FLKSTK.npflka], 
                 idisc);//<-- dummy return variable not used
        /*  Do not change these cards:*/
        Int_t igeohsm1 = 1;
        Int_t igeohsm2 = -11;
-       geohsm ( STACK.nhspnt[STACK.lstack], igeohsm1, igeohsm2, LTCLCM.mlattc );
-       STACK.nlattc[STACK.lstack] = LTCLCM.mlattc;
+       geohsm ( FLKSTK.nhspnt[FLKSTK.npflka], igeohsm1, igeohsm2, LTCLCM.mlattc );
+       FLKSTK.nlattc[FLKSTK.npflka] = LTCLCM.mlattc;
        soevsv();
     } else {
         //
index 80fb680..9183498 100644 (file)
@@ -77,13 +77,13 @@ void stupre()
     
 // Ckeck transport cut first
     Int_t ireg = EMFSTK.iremf[kp];
-    Double_t cut = (TMath::Abs(EMFSTK.ichemf[kp]) == 1) ? EFMRGN.ecut[ireg-1] :  EFMRGN.pcut[ireg-1];
+    Double_t cut = (TMath::Abs(EMFSTK.ichemf[kp]) == 1) ? EFMRGN.elethr[ireg-1] :  EFMRGN.phothr[ireg-1];
     Double_t e      = EMFSTK.etemf[kp];
     if ((e < cut) 
        && ( 
            (EMFSTK.ichemf[kp] ==  0) ||
            (EMFSTK.ichemf[kp] == -1) ||
-           (EMFSTK.ichemf[kp] ==  1 &&  EFMRGN.pcut[ireg-1] > emassmev)
+           (EMFSTK.ichemf[kp] ==  1 &&  EFMRGN.phothr[ireg-1] > emassmev)
            )
        )
     {
@@ -106,16 +106,16 @@ void stupre()
     e      *= emvgev;
     Int_t    pdg    = fluka->PDGFromId(flukaid);
     Double_t p      = sqrt(e * e - PAPROP.am[flukaid+6] * PAPROP.am[flukaid+6]);
-    Double_t px     = p * EMFSTK.u[kp];
-    Double_t pz     = p * EMFSTK.v[kp];
-    Double_t py     = p * EMFSTK.w[kp];
+    Double_t px     = p * EMFSTK.uemf[kp];
+    Double_t pz     = p * EMFSTK.vemf[kp];
+    Double_t py     = p * EMFSTK.wemf[kp];
     Double_t tof    = EMFSTK.agemf[kp];
     Double_t polx   = EMFSTK.upol[kp];
     Double_t poly   = EMFSTK.vpol[kp];
     Double_t polz   = EMFSTK.wpol[kp];
-    Double_t vx     = EMFSTK.x[kp];
-    Double_t vy     = EMFSTK.y[kp];
-    Double_t vz     = EMFSTK.z[kp];
+    Double_t vx     = EMFSTK.xemf[kp];
+    Double_t vy     = EMFSTK.yemf[kp];
+    Double_t vz     = EMFSTK.zemf[kp];
     Double_t weight = EMFSTK.wtemf[kp];
 
     Int_t ntr;
index 3562282..8eb2cdd 100644 (file)
@@ -11,9 +11,9 @@
 #include "Fdblprc.h"  //(DBLPRC) fluka common
 #include "Fevtflg.h"  //(EVTFLG) fluka common
 #include "Fpaprop.h"  //(PAPROP) fluka common
-#include "Fstack.h"   //(STACK)  fluka common
+#include "Fflkstk.h"  //(FLKSTK)  fluka common
 #include "Ftrackr.h"  //(TRACKR) fluka common
-#include "Ffinuc.h"   //(FINUC)  fluka common
+#include "Fgenstk.h"  //(GENSTK)  fluka common
 
 //Virtual MC
 #include "TFluka.h"
@@ -34,23 +34,23 @@ extern "C" {
 //*                                                                      *
 //*----------------------------------------------------------------------*
 
-// STACK.lstack  = stack pointer
-// STACK.louse   = user flag
+// FLKSTK.npflka  = stack pointer
+// FLKSTK.louse   = user flag
 // TRACKR.llouse = user defined flag for the current particle
-  STACK.louse[STACK.lstack] = TRACKR.llouse;
+  FLKSTK.louse[FLKSTK.npflka] = TRACKR.llouse;
 
 // mkbmx1 = dimension for kwb real spare array in fluka stack in DIMPAR
 // mkbmx2 = dimension for kwb int. spare array in fluka stack in DIMPAR
-// STACK.sparek  = spare real variables available for k.w.burn
-// STACK.ispark  = spare integer variables available for k.w.burn
+// FLKSTK.sparek  = spare real variables available for k.w.burn
+// FLKSTK.ispark  = spare integer variables available for k.w.burn
 // TRACKR.spausr = user defined spare variables for the current particle
 // TRACKR.ispusr = user defined spare flags for the current particle
   Int_t ispr;
   for (ispr = 0; ispr <= mkbmx1 - 1; ispr++) {
-    STACK.sparek[STACK.lstack][ispr] = TRACKR.spausr[ispr];
+    FLKSTK.sparek[FLKSTK.npflka][ispr] = TRACKR.spausr[ispr];
   }  
   for (ispr = 0; ispr <= mkbmx2 - 1; ispr++) {
-    STACK.ispark[STACK.lstack][ispr] = TRACKR.ispusr[ispr];
+    FLKSTK.ispark[FLKSTK.npflka][ispr] = TRACKR.ispusr[ispr];
   }  
  
 // Get the pointer to the VMC
@@ -71,25 +71,25 @@ extern "C" {
     Int_t done = 0;
 
     Int_t parent =  TRACKR.ispusr[mkbmx2-1];
-    Int_t kpart  = FINUC.kpart[numsec-1];
+    Int_t kpart  = GENSTK.kpart[numsec-1];
     if (kpart < -6) return;
 
     Int_t pdg = fluka->PDGFromId(kpart);
 
     
-    Double_t px = FINUC.plr[numsec-1] * FINUC.cxr[numsec-1];
-    Double_t py = FINUC.plr[numsec-1] * FINUC.cyr[numsec-1];
-    Double_t pz = FINUC.plr[numsec-1] * FINUC.czr[numsec-1];
-    Double_t e  = FINUC.tki[numsec-1] + PAPROP.am[FINUC.kpart[numsec-1]+6];
+    Double_t px = GENSTK.plr[numsec-1] * GENSTK.cxr[numsec-1];
+    Double_t py = GENSTK.plr[numsec-1] * GENSTK.cyr[numsec-1];
+    Double_t pz = GENSTK.plr[numsec-1] * GENSTK.czr[numsec-1];
+    Double_t e  = GENSTK.tki[numsec-1] + PAPROP.am[GENSTK.kpart[numsec-1]+6];
 
     Double_t vx = xx;
     Double_t vy = yy;
     Double_t vz = zz;
     
     Double_t tof  = TRACKR.atrack;
-    Double_t polx = FINUC.cxrpol[numsec-1];
-    Double_t poly = FINUC.cyrpol[numsec-1];
-    Double_t polz = FINUC.czrpol[numsec-1];
+    Double_t polx = GENSTK.cxrpol[numsec-1];
+    Double_t poly = GENSTK.cyrpol[numsec-1];
+    Double_t polz = GENSTK.czrpol[numsec-1];
     
 
     TMCProcess mech = kPHadronic;
@@ -113,7 +113,7 @@ extern "C" {
     }
     
 
-    Double_t weight = FINUC.wei[numsec-1];
+    Double_t weight = GENSTK.wei[numsec-1];
     Int_t is = 0;
     Int_t ntr;  
     // 
@@ -127,8 +127,8 @@ extern "C" {
         << numsec << "npprmr " << npprmr << endl;
 //
 //  Save current track number
-    STACK.ispark[STACK.lstack][mkbmx2-1] = ntr;
-    STACK.ispark[STACK.lstack][mkbmx2-2] = 0;
+    FLKSTK.ispark[FLKSTK.npflka][mkbmx2-1] = ntr;
+    FLKSTK.ispark[FLKSTK.npflka][mkbmx2-2] = 0;
   } // end of if (numsec-1 > npprmr)
 } // end of stuprf
 } // end of extern "C"