]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - TFluka/TFlukaConfigOption.cxx
Typo
[u/mrichter/AliRoot.git] / TFluka / TFlukaConfigOption.cxx
index fd85b313e79b35593beb15e4034219df22ebac9b..ce1360d81454e6d9fa9bd4ffaf25b344bcea7573 100644 (file)
 #include "TFlukaCerenkov.h"
 
 #include <TString.h>
+#include <TList.h>
 #include <TObjArray.h>
 #include <TVirtualMC.h>
 #include <TGeoMaterial.h>
+#include <TGeoMedium.h>
+#include <TGeoManager.h>
+#include <TGeoMedium.h>
 
 Float_t            TFlukaConfigOption::fgMatMin(-1.);
 Float_t            TFlukaConfigOption::fgMatMax(-1.);
@@ -38,11 +42,15 @@ ClassImp(TFlukaConfigOption)
 
 
 TFlukaConfigOption::TFlukaConfigOption()
+  :fMedium(-1),
+   fCMatMin(-1),
+   fCMatMax(-1),
+   fCMaterial(0)
 {
     // Default constructor
-    fMedium  = -1;
-    fCMatMin = -1.;
-    fCMatMax = -1.;    
+//    fMedium  = -1;
+//    fCMatMin = -1.;
+//    fCMatMax = -1.;    
     Int_t i;
     for (i = 0; i < 11; i++) fCutValue[i]    = -1.;
     for (i = 0; i < 15; i++) fProcessFlag[i] = -1;
@@ -50,11 +58,15 @@ TFlukaConfigOption::TFlukaConfigOption()
 
 
 TFlukaConfigOption::TFlukaConfigOption(Int_t medium)
+  :fMedium(medium),
+   fCMatMin(-1),
+   fCMatMax(-1),
+   fCMaterial(0)
 {
     // Constructor
-    fMedium = medium;
-    fCMatMin = -1.;
-    fCMatMax = -1.;    
+//    fMedium = medium;
+//    fCMatMin = -1.;
+//    fCMatMax = -1.;    
     Int_t i;
     for (i = 0; i < 11; i++) fCutValue[i]    = -1.;
     for (i = 0; i < 15; i++) fProcessFlag[i] = -1;
@@ -64,30 +76,45 @@ void TFlukaConfigOption::SetCut(const char* flagname, Double_t val)
 {
     // Set a cut value
     const TString cuts[11] = 
-       {"CUTGAM", "CUTELE", "CUTNEU", "CUTHAD", "CUTMUO", "BCUTE", "BCUTM", "DCUTE", "DCUTM", "PPCUTM", "TOFMAX"};
+       {"CUTGAM", "CUTELE", "CUTNEU", "CUTHAD", "CUTMUO", "BCUTE", "BCUTM", "DCUTE", "DCUTM", "PPCUTM", "TOFMAX"};
     Int_t i;
     for (i = 0; i < 11; i++) {
-       if (cuts[i].CompareTo(flagname) == 0) {
-           fCutValue[i] = val;
-           if (fMedium == -1) fgDCutValue[i] = val;
-           break;
-       }
+       if (cuts[i].CompareTo(flagname) == 0) {
+           fCutValue[i] = val;
+           if (fMedium == -1) fgDCutValue[i] = val;
+           break;
+       }
     }
 }
 
+void TFlukaConfigOption::SetModelParameter(const char* flagname, Double_t val)
+{
+    // Set a model parameter value
+    const TString parms[2] = {"PRIMIO_N", "PRIMIO_E"};
+    Int_t i;
+    for (i = 0; i < 2; i++) {
+       if (parms[i].CompareTo(flagname) == 0) {
+           fModelParameter[i] = val;
+           break;
+       }
+    }
+}
+
+
 void TFlukaConfigOption::SetProcess(const char* flagname, Int_t flag)
 {
     // Set a process flag
     const TString process[15] = 
-       {"DCAY", "PAIR", "COMP", "PHOT", "PFIS", "DRAY", "ANNI", "BREM", "MUNU", "CKOV", 
-        "HADR", "LOSS", "MULS", "RAYL", "STRA"};
+       {"DCAY", "PAIR", "COMP", "PHOT", "PFIS", "DRAY", "ANNI", "BREM", "MUNU", "CKOV",
+        "HADR", "LOSS", "MULS", "RAYL", "STRA"};
+    
     Int_t i;
     for (i = 0; i < 15; i++) {
-       if (process[i].CompareTo(flagname) == 0) {
-           fProcessFlag[i] = flag;
-           if (fMedium == -1) fgDProcessFlag[i] = flag;
-           break;
-       }
+       if (process[i].CompareTo(flagname) == 0) {
+           fProcessFlag[i] = flag;
+           if (fMedium == -1) fgDProcessFlag[i] = flag;
+           break;
+       }
     }
 }
 
@@ -96,19 +123,75 @@ void TFlukaConfigOption::WriteFlukaInputCards()
     // Write the FLUKA input cards for the set of process flags and cuts
     //
     //
-    if (fMedium > -1) {
-       fprintf(fgFile,"*\n*Material specific process and cut settings for #%8d \n", fMedium);
-       fCMatMin = fMedium;
-       fCMatMax = fMedium;
+    //
+    // Check if global option or medium specific
+
+    Bool_t mediumIsSensitive = kFALSE;
+    TGeoMedium*   med    = 0x0;
+    TGeoMedium*   medium = 0x0;    
+    TGeoMaterial* mat    = 0x0;
+
+    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) {
+           // Material not used -- return
+           return;
+       }
+       //
+       // Find material
+       TObjArray *matList = fluka->GetFlukaMaterials();
+       Int_t nmaterial =  matList->GetEntriesFast();
+       fCMaterial = 0;
+       for (Int_t im = 0; im < nmaterial; im++)
+       {
+           fCMaterial = dynamic_cast<TGeoMaterial*> (matList->At(im));
+           Int_t idmat = fCMaterial->GetIndex();
+           if (idmat == fMedium) break;
+       } // materials
+        //
+       // Find medium
+       TList *medlist = gGeoManager->GetListOfMedia();
+       TIter next(medlist);
+       while((med = (TGeoMedium*)next()))
+       {
+           mat = med->GetMaterial();
+           if (mat->GetIndex() == fMedium) {
+              medium = med;
+              break;
+           }
+       } // media
+       //
+       // Check if sensitive
+       if (medium->GetParam(0) != 0.) mediumIsSensitive = kTRUE;
+
+       fprintf(fgFile,"*\n*Material specific process and cut settings for #%8d %s\n", fMedium, fCMaterial->GetName());
+       fCMatMin = fMedium;
+       fCMatMax = fMedium;
     } else {
-       fprintf(fgFile,"*\n*Global process and cut settings \n");
-       fCMatMin = fgMatMin;
-       fCMatMax = fgMatMax;
+       fprintf(fgFile,"*\n*Global process and cut settings \n");
+       fCMatMin = fgMatMin;
+       fCMatMax = fgMatMax;
     }
 
 //
 // Handle Process Flags 
+//
+//
+//  First make sure that all cuts are taken into account
+    if (DefaultProcessFlag(kPAIR) > 0 && fProcessFlag[kPAIR] == -1 && (fCutValue[kCUTELE] >= 0. || fCutValue[kPPCUTM] >= 0.)) 
+       fProcessFlag[kPAIR] = DefaultProcessFlag(kPAIR);
+    if (DefaultProcessFlag(kBREM) > 0 && fProcessFlag[kBREM] == -1 && (fCutValue[kBCUTE]  >= 0. || fCutValue[kBCUTM] >= 0.)) 
+       fProcessFlag[kBREM] = DefaultProcessFlag(kBREM);
+    if (DefaultProcessFlag(kDRAY) > 0 && fProcessFlag[kDRAY] == -1 && (fCutValue[kDCUTE]  >= 0. || fCutValue[kDCUTM] >= 0.)) 
+       fProcessFlag[kDRAY] = DefaultProcessFlag(kDRAY);
 //    
+//
     if (fProcessFlag[kDCAY] != -1) ProcessDCAY();
     if (fProcessFlag[kPAIR] != -1) ProcessPAIR();
     if (fProcessFlag[kBREM] != -1) ProcessBREM();
@@ -132,8 +215,13 @@ void TFlukaConfigOption::WriteFlukaInputCards()
     if (fCutValue[kCUTNEU] >= 0.) ProcessCUTNEU();
     if (fCutValue[kCUTHAD] >= 0.) ProcessCUTHAD();
     if (fCutValue[kCUTMUO] >= 0.) ProcessCUTMUO();
-
+//
+//  Time of flight 
     if (fCutValue[kTOFMAX] >= 0.) ProcessTOFMAX();
+
+//
+//  Tracking precission
+    if (mediumIsSensitive) ProcessSensitiveMedium();
 }
 
 void TFlukaConfigOption::ProcessDCAY()
@@ -141,9 +229,9 @@ void TFlukaConfigOption::ProcessDCAY()
     // Process DCAY option
     fprintf(fgFile,"*\n* --- DCAY --- Decays. Flag = %5d\n", fProcessFlag[kDCAY]);
     if (fProcessFlag[kDCAY] == 0) {
-       printf("Decays cannot be switched off \n");
+       printf("Decays cannot be switched off \n");
     } else {
-       fprintf(fgFile, "* Decays are on by default\n");
+       fprintf(fgFile, "* Decays are on by default\n");
     }
 }
 
@@ -151,15 +239,15 @@ 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.);
+       fprintf(fgFile,"EMFCUT    %10.1f%10.1f%10.4g%10.1f%10.1f%10.1fPHOT-THR\n",0., 0., 1e10,  fCMatMin, fCMatMax, 1.);
     }
     
     //
@@ -170,7 +258,7 @@ void TFlukaConfigOption::ProcessPAIR()
     //
 
     if (fProcessFlag[kBREM] == -1 ) fProcessFlag[kBREM] = fgDProcessFlag[kBREM];
-    if (fCutValue[kBCUTM]   == -1.) fCutValue[kBCUTM]   = fgDCutValue[kBCUTM]; 
+    if (fCutValue[kBCUTM]   == -1.) fCutValue[kBCUTM]   = fgDCutValue[kBCUTM];       
 
 
     Float_t flag = -3.;    
@@ -179,13 +267,13 @@ void TFlukaConfigOption::ProcessPAIR()
     if (fProcessFlag[kPAIR] >  0 && fProcessFlag[kBREM]  > 0) flag =  3.;    
     if (fProcessFlag[kPAIR] == 0 && fProcessFlag[kBREM] == 0) flag = -3.;
     // Flag BREM card as handled
-    fProcessFlag[kBREM] = -1;
+    fProcessFlag[kBREM] = - fProcessFlag[kBREM];
     
     //
     // Energy cut for pair prodution
     //
     Float_t cutP = fCutValue[kPPCUTM];
-    if (fCutValue[kPPCUTM]   == -1.) cutP = fgDCutValue[kPPCUTM];      
+    if (fCutValue[kPPCUTM]   == -1.) cutP = fgDCutValue[kPPCUTM];       
     // In G3 this is the cut on the total energy of the e+e- pair
     // In FLUKA the cut is on the kinetic energy of the electron and poistron
     cutP = cutP / 2. - 0.51099906e-3;
@@ -197,12 +285,12 @@ void TFlukaConfigOption::ProcessPAIR()
     //
     Float_t cutB = 0.;
     if (flag > 1.) {
-       fprintf(fgFile,"*\n* +++ BREM --- Bremsstrahlung by muons/hadrons. Flag = %5d, BCUTM = %13.4g \n",
-           fProcessFlag[kBREM], fCutValue[kBCUTM]);
+       fprintf(fgFile,"*\n* +++ BREM --- Bremsstrahlung by muons/hadrons. Flag = %5d, BCUTM = %13.4g \n",
+           fProcessFlag[kBREM], fCutValue[kBCUTM]);
 
-       cutB =  fCutValue[kBCUTM];
-       // No explicite production of gammas
-       if (fProcessFlag[kBREM] == 2) cutB = -1.;
+       cutB =  fCutValue[kBCUTM];
+       // No explicite production of gammas
+       if (fProcessFlag[kBREM] == 2) cutB = -1.;
     }
 
     fprintf(fgFile,"PAIRBREM  %10.1f%10.4g%10.4g%10.1f%10.1f\n",flag, cutP, cutB, fCMatMin, fCMatMax);
@@ -213,26 +301,26 @@ void TFlukaConfigOption::ProcessBREM()
 {
     // Process BREM option
     fprintf(fgFile,"*\n* --- BREM --- Bremsstrahlung by e+/- and muons/hadrons. Flag = %5d, BCUTE = %13.4g, BCUTM = %13.4g \n",
-           fProcessFlag[kBREM], fCutValue[kBCUTE], fCutValue[kBCUTM]);
+           fProcessFlag[kBREM], fCutValue[kBCUTE], fCutValue[kBCUTM]);
 
     //
     // e+/- -> e+/- gamma
     //
     Float_t cutB = fCutValue[kBCUTE];
-    if (fCutValue[kBCUTE]   == -1.) cutB = fgDCutValue[kBCUTE];        
+    if (fCutValue[kBCUTE]   == -1.) cutB = fgDCutValue[kBCUTE];       
     
     
-    if (fProcessFlag[kBREM] > 0) {
-       fprintf(fgFile,"EMFCUT    %10.4g%10.1f%10.1f%10.1f%10.1f%10.1fELPO-THR\n",cutB,  0., 0.,  fCMatMin, fCMatMax, 1.);
+    if (TMath::Abs(fProcessFlag[kBREM]) > 0) {
+       fprintf(fgFile,"EMFCUT    %10.4g%10.1f%10.1f%10.1f%10.1f%10.1fELPO-THR\n",cutB,  0., 0.,  fCMatMin, fCMatMax, 1.);
     } else {
-       fprintf(fgFile,"EMFCUT    %10.4g%10.1f%10.1f%10.1f%10.1f%10.1fELPO-THR\n",1.e10, 0., 0.,  fCMatMin, fCMatMax, 1.);
+       fprintf(fgFile,"EMFCUT    %10.4g%10.1f%10.1f%10.1f%10.1f%10.1fELPO-THR\n",1.e10, 0., 0.,  fCMatMin, fCMatMax, 1.);
     }
 
     //
     // Bremsstrahlung by muons and hadrons
     //
     cutB = fCutValue[kBCUTM];
-    if (fCutValue[kBCUTM]   == -1.) cutB = fgDCutValue[kBCUTM];        
+    if (fCutValue[kBCUTM]   == -1.) cutB = fgDCutValue[kBCUTM];       
     if (fProcessFlag[kBREM] == 2) cutB = -1.;
     Float_t flag = 2.;
     if (fProcessFlag[kBREM] == 0) flag = -2.;
@@ -250,9 +338,9 @@ void TFlukaConfigOption::ProcessCOMP()
     //
 
     if (fProcessFlag[kCOMP] > 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.1f%10.1f%10.1f%10.1fPHOT-THR\n",0.   , 0., 0.,  fCMatMin, fCMatMax, 1.);
     } else {
-       fprintf(fgFile,"EMFCUT    %10.4g%10.1f%10.1f%10.1f%10.1f%10.1fPHOT-THR\n",1.e10, 0., 0.,  fCMatMin, fCMatMax, 1.);
+       fprintf(fgFile,"EMFCUT    %10.4g%10.1f%10.1f%10.1f%10.1f%10.1fPHOT-THR\n",1.e10, 0., 0.,  fCMatMin, fCMatMax, 1.);
     }
 }
 
@@ -266,9 +354,9 @@ 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.);
+       fprintf(fgFile,"EMFCUT    %10.1f%10.4g%10.1f%10.1f%10.1f%10.1fPHOT-THR\n",0., 1.e10, 0.,  fCMatMin, fCMatMax, 1.);
     }
 }
 
@@ -282,9 +370,9 @@ void TFlukaConfigOption::ProcessANNI()
     //
 
     if (fProcessFlag[kANNI] > 0) {
-       fprintf(fgFile,"EMFCUT    %10.1f%10.1f%10.1f%10.1f%10.1f%10.1fANNH-THR\n",0.   , 0., 0.,  fCMatMin, fCMatMax, 1.);
+       fprintf(fgFile,"EMFCUT    %10.1f%10.1f%10.1f%10.1f%10.1f%10.1fANNH-THR\n",0.   , 0., 0.,  fCMatMin, fCMatMax, 1.);
     } else {
-       fprintf(fgFile,"EMFCUT    %10.4g%10.1f%10.1f%10.1f%10.1f%10.1fANNH-THR\n",1.e10, 0., 0.,  fCMatMin, fCMatMax, 1.);
+       fprintf(fgFile,"EMFCUT    %10.4g%10.1f%10.1f%10.1f%10.1f%10.1fANNH-THR\n",1.e10, 0., 0.,  fCMatMin, fCMatMax, 1.);
     }
 }
 
@@ -299,9 +387,9 @@ void TFlukaConfigOption::ProcessPFIS()
     //
 
     if (fProcessFlag[kPFIS] > 0) {
-       fprintf(fgFile,"PHOTONUC  %10.1f%10.1f%10.1f%10.1f%10.1f%10.1f\n",(Float_t) fProcessFlag[kPFIS], 0., 0.,  fCMatMin, fCMatMax, 1.);
+       fprintf(fgFile,"PHOTONUC  %10.1f%10.1f%10.1f%10.1f%10.1f%10.1f\n",(Float_t) fProcessFlag[kPFIS], 0., 0.,  fCMatMin, fCMatMax, 1.);
     } else {
-       fprintf(fgFile,"PHOTONUC  %10.1f%10.1f%10.1f%10.1f%10.1f%10.1f\n",-1.                          , 0., 0.,  fCMatMin, fCMatMax, 1.);
+       fprintf(fgFile,"PHOTONUC  %10.1f%10.1f%10.1f%10.1f%10.1f%10.1f\n",-1.                          , 0., 0.,  fCMatMin, fCMatMax, 1.);
     }
 }
 
@@ -314,9 +402,9 @@ void TFlukaConfigOption::ProcessMUNU()
     // Muon nuclear interactions
     //
     if (fProcessFlag[kMUNU] > 0) {
-       fprintf(fgFile,"MUPHOTON  %10.1f%10.3f%10.3f%10.1f%10.1f%10.1f\n",(Float_t )fProcessFlag[kMUNU], 0.25, 0.75,  fCMatMin, fCMatMax, 1.);
+       fprintf(fgFile,"MUPHOTON  %10.1f%10.3f%10.3f%10.1f%10.1f%10.1f\n",(Float_t )fProcessFlag[kMUNU], 0.25, 0.75,  fCMatMin, fCMatMax, 1.);
     } else {
-       fprintf(fgFile,"MUPHOTON  %10.1f%10.1f%10.1f%10.1f%10.1f%10.1f\n",-1.                          , 0.,   0.,    fCMatMin, fCMatMax, 1.);
+       fprintf(fgFile,"MUPHOTON  %10.1f%10.1f%10.1f%10.1f%10.1f%10.1f\n",-1.                          , 0.,   0.,    fCMatMin, fCMatMax, 1.);
     }
 }
 
@@ -333,12 +421,12 @@ void TFlukaConfigOption::ProcessRAYL()
     //
     // Loop over regions of a given material
     for (Int_t k = 0; k < nreg; k++) {
-       Float_t ireg = reglist[k];
-       if (fProcessFlag[kRAYL] > 0) {
-           fprintf(fgFile,"EMFRAY    %10.1f%10.1f%10.1f%10.1f\n", 1., ireg, ireg, 1.);
-       } else {
-           fprintf(fgFile,"EMFRAY    %10.1f%10.1f%10.1f%10.1f\n", 3., ireg, ireg, 1.);
-       }
+       Float_t ireg = reglist[k];
+       if (fProcessFlag[kRAYL] > 0) {
+           fprintf(fgFile,"EMFRAY    %10.1f%10.1f%10.1f%10.1f\n", 1., ireg, ireg, 1.);
+       } else {
+           fprintf(fgFile,"EMFRAY    %10.1f%10.1f%10.1f%10.1f\n", 3., ireg, ireg, 1.);
+       }
     }
 }
 
@@ -356,54 +444,54 @@ void TFlukaConfigOption::ProcessCKOV()
     Int_t nmaterial =  matList->GetEntriesFast();
     for (Int_t im = 0; im < nmaterial; im++)
     {
-       TGeoMaterial* material = dynamic_cast<TGeoMaterial*> (matList->At(im));
-       Int_t idmat = material->GetIndex();
+       TGeoMaterial* material = dynamic_cast<TGeoMaterial*> (matList->At(im));
+       Int_t idmat = material->GetIndex();
 //
 // Check if global option
-       if (fMedium != -1 && idmat != fMedium) continue;
-       
-       TFlukaCerenkov* cerenkovProp;
-       if (!(cerenkovProp = dynamic_cast<TFlukaCerenkov*>(material->GetCerenkovProperties()))) continue;
-       //
-       // This medium has Cerenkov properties 
-       //
-       //
-       if (fMedium == -1 || (fMedium != -1 && fProcessFlag[kCKOV] > 0)) {
-           // Write OPT-PROD card for each medium 
-           Float_t  emin  = cerenkovProp->GetMinimumEnergy();
-           Float_t  emax  = cerenkovProp->GetMaximumEnergy();        
-           fprintf(fgFile, "OPT-PROD  %10.4g%10.4g%10.4g%10.4g%10.4g%10.4gCERENKOV\n", emin, emax, 0., 
-                   Float_t(idmat), Float_t(idmat), 0.); 
-           //
-           // Write OPT-PROP card for each medium 
-           // Forcing FLUKA to call user routines (queffc.cxx, rflctv.cxx, rfrndx.cxx)
-           //
-           fprintf(fgFile, "OPT-PROP  %10.4g%10.4g%10.4g%10.1f%10.1f%10.1fWV-LIMIT\n",  
-                   cerenkovProp->GetMinimumWavelength(), cerenkovProp->GetMaximumWavelength(), cerenkovProp->GetMaximumWavelength(), 
-                   Float_t(idmat), Float_t(idmat), 0.0);
-           
-
-           fprintf(fgFile, "OPT-PROP  %10.1f%10.1f%10.1f%10.1f%10.1f%10.1f\n", -100., -100., -100., 
-                   Float_t(idmat), Float_t(idmat), 0.0);
-           
-           for (Int_t j = 0; j < 3; j++) {
-               fprintf(fgFile, "OPT-PROP  %10.1f%10.1f%10.1f%10.1f%10.1f%10.1f&\n", -100., -100., -100., 
-                       Float_t(idmat), Float_t(idmat), 0.0);
-           }
-
-
-           // Photon detection efficiency user defined     
-           if (cerenkovProp->IsSensitive())
-               fprintf(fgFile, "OPT-PROP  %10.1f%10.1f%10.1f%10.1f%10.1f%10.1fSENSITIV\n", -100., -100., -100., 
-                       Float_t(idmat), Float_t(idmat), 0.0);
-           // Material has a reflective surface
-           if (cerenkovProp->IsMetal())
-               fprintf(fgFile, "OPT-PROP  %10.1f%10.1f%10.1f%10.1f%10.1f%10.1fMETAL\n", -100., -100., -100., 
-                       Float_t(idmat), Float_t(idmat), 0.0);
-
-       } else {
-           fprintf(fgFile,"OPT-PROD  %10.1f%10.1f%10.1f%10.1f%10.1f%10.1fCERE-OFF\n",0., 0., 0., fCMatMin, fCMatMax, 1.);
-       }
+       if (fMedium != -1 && idmat != fMedium) continue;
+       
+       TFlukaCerenkov* cerenkovProp;
+       if (!(cerenkovProp = dynamic_cast<TFlukaCerenkov*>(material->GetCerenkovProperties()))) continue;
+       //
+       // This medium has Cerenkov properties
+       //
+       //
+       if (fMedium == -1 || (fMedium != -1 && fProcessFlag[kCKOV] > 0)) {
+           // Write OPT-PROD card for each medium
+           Float_t  emin  = cerenkovProp->GetMinimumEnergy();
+           Float_t  emax  = cerenkovProp->GetMaximumEnergy();
+           fprintf(fgFile, "OPT-PROD  %10.4g%10.4g%10.4g%10.4g%10.4g%10.4gCERENKOV\n", emin, emax, 0.,
+                  Float_t(idmat), Float_t(idmat), 0.);
+           //
+           // Write OPT-PROP card for each medium
+           // Forcing FLUKA to call user routines (queffc.cxx, rflctv.cxx, rfrndx.cxx)
+           //
+           fprintf(fgFile, "OPT-PROP  %10.4g%10.4g%10.4g%10.1f%10.1f%10.1fWV-LIMIT\n",
+                  cerenkovProp->GetMinimumWavelength(), cerenkovProp->GetMaximumWavelength(), cerenkovProp->GetMaximumWavelength(),
+                  Float_t(idmat), Float_t(idmat), 0.0);
+       
+
+           fprintf(fgFile, "OPT-PROP  %10.1f%10.1f%10.1f%10.1f%10.1f%10.1f\n", -100., -100., -100.,
+                  Float_t(idmat), Float_t(idmat), 0.0);
+       
+           for (Int_t j = 0; j < 3; j++) {
+              fprintf(fgFile, "OPT-PROP  %10.1f%10.1f%10.1f%10.1f%10.1f%10.1f&\n", -100., -100., -100.,
+                     Float_t(idmat), Float_t(idmat), 0.0);
+           }
+
+
+           // Photon detection efficiency user defined
+           if (cerenkovProp->IsSensitive())
+              fprintf(fgFile, "OPT-PROP  %10.1f%10.1f%10.1f%10.1f%10.1f%10.1fSENSITIV\n", -100., -100., -100.,
+                     Float_t(idmat), Float_t(idmat), 0.0);
+           // Material has a reflective surface
+           if (cerenkovProp->IsMetal())
+              fprintf(fgFile, "OPT-PROP  %10.1f%10.1f%10.1f%10.1f%10.1f%10.1fMETAL\n", -100., -100., -100.,
+                     Float_t(idmat), Float_t(idmat), 0.0);
+
+       } else {
+           fprintf(fgFile,"OPT-PROD  %10.1f%10.1f%10.1f%10.1f%10.1f%10.1fCERE-OFF\n",0., 0., 0., fCMatMin, fCMatMax, 1.);
+       }
     }
 }
 
@@ -414,10 +502,13 @@ void TFlukaConfigOption::ProcessHADR()
     fprintf(fgFile,"*\n* --- HADR --- Hadronic interactions. Flag = %5d\n", fProcessFlag[kHADR]);
     
     if (fProcessFlag[kHADR] > 0) {
-       fprintf(fgFile,"*\n*Hadronic interaction is ON by default in FLUKA\n");
+       fprintf(fgFile,"*\n*Hadronic interaction is ON by default in FLUKA\n");
     } else {
-       if (fMedium != -1) printf("Hadronic interactions cannot be switched off material by material !\n");
-       fprintf(fgFile,"THRESHOL  %10.1f%10.1f%10.1f%10.1e%10.1f\n",0., 0., 0., 1.e10, 0.);
+       if (fMedium != -1) {
+          printf("Hadronic interactions cannot be switched off material by material !\n");
+       } else {
+          fprintf(fgFile,"THRESHOL  %10.1f%10.1f%10.1f%10.1e%10.1f\n",0., 0., 0., 1.e10, 0.);
+       }
     }
 }
 
@@ -431,9 +522,9 @@ void TFlukaConfigOption::ProcessMULS()
     // Multiple scattering
     //
     if (fProcessFlag[kMULS] > 0) {
-       fprintf(fgFile,"*\n*Multiple scattering is  ON by default in FLUKA\n");
+       fprintf(fgFile,"*\n*Multiple scattering is  ON by default in FLUKA\n");
     } else {
-       fprintf(fgFile,"MULSOPT   %10.1f%10.1f%10.1f%10.1f%10.1f\n",0., 3., 3., fCMatMin, fCMatMax);
+       fprintf(fgFile,"MULSOPT   %10.1f%10.1f%10.1f%10.1f%10.1f\n",0., 3., 3., fCMatMin, fCMatMax);
     }
 }
 
@@ -441,34 +532,36 @@ void TFlukaConfigOption::ProcessLOSS()
 {
     // Process LOSS option
     fprintf(fgFile,"*\n* --- LOSS --- Ionisation energy loss. Flags: LOSS = %5d, DRAY = %5d, STRA = %5d; Cuts: DCUTE = %13.4g, DCUTM = %13.4g \n",
-           fProcessFlag[kLOSS], fProcessFlag[kDRAY], fProcessFlag[kSTRA], fCutValue[kDCUTE], fCutValue[kDCUTM]);
+           fProcessFlag[kLOSS], fProcessFlag[kDRAY], fProcessFlag[kSTRA], fCutValue[kDCUTE], fCutValue[kDCUTM]);
     //
     // Ionisation energy loss
     //
     //
     // Impose consistency
     
-    if (fProcessFlag[kLOSS] == 1 || fProcessFlag[kLOSS] == 3) {
-       fProcessFlag[kDRAY] = 1;
+    if (fProcessFlag[kLOSS] == 1 || fProcessFlag[kLOSS] == 3 || fProcessFlag[kLOSS] > 10) {
+    // Restricted fluctuations
+       fProcessFlag[kDRAY] = 1;
     } else if (fProcessFlag[kLOSS] == 2) {
-       fProcessFlag[kDRAY] = 0;
-       fCutValue[kDCUTE] = 1.e10;
-       fCutValue[kDCUTM] = 1.e10;      
+    // Full fluctuations
+       fProcessFlag[kDRAY] = 0;
+       fCutValue[kDCUTE] = 1.e10;
+       fCutValue[kDCUTM] = 1.e10;
     } else {
-       if (fProcessFlag[kDRAY] == 1) {
-           fProcessFlag[kLOSS] = 1;
-       } else if (fProcessFlag[kDRAY] == 0) {
-           fProcessFlag[kLOSS] = 2;
-           fCutValue[kDCUTE] = 1.e10;
-           fCutValue[kDCUTM] = 1.e10;  
-       }
+       if (fProcessFlag[kDRAY] == 1) {
+           fProcessFlag[kLOSS] = 1;
+       } else if (fProcessFlag[kDRAY] == 0) {
+           fProcessFlag[kLOSS] = 2;
+           fCutValue[kDCUTE] = 1.e10;
+           fCutValue[kDCUTM] = 1.e10;
+       }
     }
     
     if (fCutValue[kDCUTE] == -1.) fCutValue[kDCUTE] = fgDCutValue[kDCUTE];
     if (fCutValue[kDCUTM] == -1.) fCutValue[kDCUTM] = fgDCutValue[kDCUTM];    
     
     Float_t cutM =  fCutValue[kDCUTM];    
-       
+       
 
     if (fProcessFlag[kSTRA] == -1) fProcessFlag[kSTRA] = fgDProcessFlag[kSTRA];
     if (fProcessFlag[kSTRA] ==  0) fProcessFlag[kSTRA] = 1;
@@ -479,20 +572,34 @@ void TFlukaConfigOption::ProcessLOSS()
 //
 // Restricted energy loss fluctuations 
 //
-       fprintf(fgFile,"IONFLUCT  %10.1f%10.1f%10.1f%10.1f%10.1f\n", 1., 1., stra, fCMatMin, fCMatMax);
-       fprintf(fgFile,"DELTARAY  %10.4g%10.1f%10.1f%10.1f%10.1f%10.1f\n", cutM, 0., 0., fCMatMin, fCMatMax, 1.);
+       fprintf(fgFile,"IONFLUCT  %10.1f%10.1f%10.1f%10.1f%10.1f\n", 1., 1., stra, fCMatMin, fCMatMax);
+       fprintf(fgFile,"DELTARAY  %10.4g%10.1f%10.1f%10.1f%10.1f%10.1f\n", cutM, 0., 0., fCMatMin, fCMatMax, 1.);
+    } else if (fProcessFlag[kLOSS] > 10) {
+//
+// Primary ionisation electron generation
+//
+       // Ionisation model
+       Float_t ioModel = Float_t (fProcessFlag[kLOSS]-10);
+       //  Effective 1st ionisation potential
+       Float_t ePot    = ModelParameter(kPRIMIOE);
+       // Number of primary ionisations per cm for a mip
+       Float_t nPrim   = ModelParameter(kPRIMION);
+       
+       fprintf(fgFile,"IONFLUCT  %10.1f%10.1f%10.1f%10.1f%10.1f\n", 1., 1., stra, fCMatMin, fCMatMax);
+       fprintf(fgFile,"IONFLUCT  %10.1f%10.1f%10.1f%10.1f%10.1f%10.1fPRIM-ION\n", ePot, nPrim, ioModel, fCMatMin, fCMatMax, 1.);
+       fprintf(fgFile,"DELTARAY  %10.4g%10.1f%10.1f%10.1f%10.1f%10.1f\n", cutM, 0., 0., fCMatMin, fCMatMax, 1.);
     } else if (fProcessFlag[kLOSS] == 4) {
 //
 // No fluctuations
 //
-       fprintf(fgFile,"IONFLUCT  %10.1f%10.1f%10.1f%10.1f%10.1f\n",-1., -1., stra, fCMatMin, fCMatMax);        
-       fprintf(fgFile,"DELTARAY  %10.4g%10.1f%10.1f%10.1f%10.1f%10.1f\n", 1.e10, 0., 0., fCMatMin, fCMatMax, 1.);      
+       fprintf(fgFile,"IONFLUCT  %10.1f%10.1f%10.1f%10.1f%10.1f\n",-1., -1., stra, fCMatMin, fCMatMax);
+       fprintf(fgFile,"DELTARAY  %10.4g%10.1f%10.1f%10.1f%10.1f%10.1f\n", 1.e10, 0., 0., fCMatMin, fCMatMax, 1.);
     }  else {
 //
 // Full fluctuations
 //
-       fprintf(fgFile,"IONFLUCT  %10.1f%10.1f%10.1f%10.1f%10.1f\n",1., 1., stra, fCMatMin, fCMatMax);  
-       fprintf(fgFile,"DELTARAY  %10.4g%10.1f%10.1f%10.1f%10.1f%10.1f\n", 1.e10, 0., 0., fCMatMin, fCMatMax, 1.);      
+       fprintf(fgFile,"IONFLUCT  %10.1f%10.1f%10.1f%10.1f%10.1f\n",1., 1., stra, fCMatMin, fCMatMax);
+       fprintf(fgFile,"DELTARAY  %10.4g%10.1f%10.1f%10.1f%10.1f%10.1f\n", 1.e10, 0., 0., fCMatMin, fCMatMax, 1.);
     }
 }
 
@@ -503,18 +610,27 @@ void TFlukaConfigOption::ProcessCUTGAM()
 //
     fprintf(fgFile,"*\n*Cut for Gammas. CUTGAM = %13.4g\n", fCutValue[kCUTGAM]);
     if (fMedium == -1) {
-       fprintf(fgFile,"EMFCUT    %10.4g%10.4g%10.1f%10.1f%10.1f%10.1f\n", 
-               0., fCutValue[kCUTGAM], 0., 0., Float_t(fgGeom->NofVolumes()), 1.);
+       fprintf(fgFile,"EMFCUT    %10.4g%10.4g%10.1f%10.1f%10.1f%10.1f\n",
+              0., fCutValue[kCUTGAM], 0., 0., Float_t(fgGeom->NofVolumes()), 1.);
+
     } else {
-       Int_t nreg, *reglist;
-       Float_t ireg;
-       reglist = fgGeom->GetMaterialList(fMedium, nreg);
-       // Loop over regions of a given material
-       for (Int_t k = 0; k < nreg; k++) {
-           ireg = reglist[k];
-           fprintf(fgFile,"EMFCUT    %10.4g%10.4g%10.1f%10.1f%10.1f%10.1f\n", 0.,fCutValue[kCUTGAM], 0., ireg, ireg, 1.);
-       }
+       Int_t nreg, *reglist;
+       Float_t ireg;
+       reglist = fgGeom->GetMaterialList(fMedium, nreg);
+       // Loop over regions of a given material
+       for (Int_t k = 0; k < nreg; k++) {
+           ireg = reglist[k];
+           fprintf(fgFile,"EMFCUT    %10.4g%10.4g%10.1f%10.1f%10.1f%10.1f\n", 0.,fCutValue[kCUTGAM], 0., ireg, ireg, 1.);
+       }
     }
+    
+    // Transport production cut used for pemf
+    //
+    //  FUDGEM paramter. The parameter takes into account th contribution of atomic electrons to multiple scattering.
+    //  For production and transport cut-offs larger than 100 keV it must be set = 1.0, while in the keV region it must be
+    Float_t parFudgem = (fCutValue[kCUTGAM] > 1.e-4)? 1.0 : 0.0 ;
+    fprintf(fgFile,"EMFCUT    %10.4g%10.4g%10.1f%10.1f%10.1f%10.1fPROD-CUT\n", 
+           0., fCutValue[kCUTGAM], parFudgem, fCMatMin, fCMatMax, 1.);
 }
 
 void TFlukaConfigOption::ProcessCUTELE()
@@ -523,69 +639,115 @@ void TFlukaConfigOption::ProcessCUTELE()
 //
     fprintf(fgFile,"*\n*Cut for e+/e-. CUTELE = %13.4g\n", fCutValue[kCUTELE]);
     if (fMedium == -1) {
-       fprintf(fgFile,"EMFCUT    %10.4g%10.4g%10.1f%10.1f%10.1f%10.1f\n", 
-               -fCutValue[kCUTELE], 0., 0., 0., Float_t(fgGeom->NofVolumes()), 1.);
+       fprintf(fgFile,"EMFCUT    %10.4g%10.4g%10.1f%10.1f%10.1f%10.1f\n",
+              -fCutValue[kCUTELE], 0., 0., 0., Float_t(fgGeom->NofVolumes()), 1.);
     } else {
-       Int_t nreg, *reglist;
-       Float_t ireg;
-       reglist = fgGeom->GetMaterialList(fMedium, nreg);
-       // Loop over regions of a given material
-       for (Int_t k = 0; k < nreg; k++) {
-           ireg = reglist[k];
-           fprintf(fgFile,"EMFCUT    %10.4g%10.4g%10.1f%10.1f%10.1f%10.1f\n", -fCutValue[kCUTELE], 0., 0., ireg, ireg, 1.);
-       }
+       Int_t nreg, *reglist;
+       Float_t ireg;
+       reglist = fgGeom->GetMaterialList(fMedium, nreg);
+       // Loop over regions of a given material
+       for (Int_t k = 0; k < nreg; k++) {
+           ireg = reglist[k];
+           fprintf(fgFile,"EMFCUT    %10.4g%10.4g%10.1f%10.1f%10.1f%10.1f\n", -fCutValue[kCUTELE], 0., 0., ireg, ireg, 1.);
+       }
     }
+    // Transport production cut used for pemf
+    //
+    //  FUDGEM paramter. The parameter takes into account th contribution of atomic electrons to multiple scattering.
+    //  For production and transport cut-offs larger than 100 keV it must be set = 1.0, while in the keV region it must be
+    Float_t parFudgem = (fCutValue[kCUTELE] > 1.e-4)? 1.0 : 0.0;
+    fprintf(fgFile,"EMFCUT    %10.4g%10.4g%10.1f%10.1f%10.1f%10.1fPROD-CUT\n", 
+           -fCutValue[kCUTELE], 0., parFudgem, fCMatMin, fCMatMax, 1.);
 }
 
 void TFlukaConfigOption::ProcessCUTNEU()
 {
-    // Cut on neutral hadrons
-    fprintf(fgFile,"*\n*Cut for neutral hadrons. CUTNEU = %13.4g\n", fCutValue[kCUTNEU]);
-    if (fMedium == -1) {
-       Float_t cut = fCutValue[kCUTNEU];
-       //
-       // 8.0 = Neutron
-       // 9.0 = Antineutron
-       //
-       // If the cut is > 19.6 MeV it is assumed the low energy neutron transport is requested.
-       // In this case the cut has to coincide with the upper  limit of the first energy group.
-       //
-       Float_t neutronCut = cut;
-       if (neutronCut < 0.0196) {
-           neutronCut = 0.0196;
-           printf("Cut on neutron lower than upper limit if first energy group.\n");
-           printf("Cut reset to 19.6 MeV !\n");
-       }
-       fprintf(fgFile,"PART-THR  %10.4g%10.1f%10.1f\n", -neutronCut,  8.0,  9.0);
-       //
-       // 12.0 = Kaon zero long
-       fprintf(fgFile,"PART-THR  %10.4g%10.1f%10.1f\n", -cut, 12.0, 12.0);
-       // 17.0 = Lambda, 18.0 = Antilambda
-       // 19.0 = Kaon zero short
-       fprintf(fgFile,"PART-THR  %10.4g%10.1f%10.1f\n", -cut, 17.0, 19.0);
-       // 22.0 = Sigma zero, Pion zero, Kaon zero
-       // 25.0 = Antikaon zero
-       fprintf(fgFile,"PART-THR  %10.4g%10.1f%10.1f\n", -cut, 22.0, 25.0);
-       // 32.0 = Antisigma zero
-       fprintf(fgFile,"PART-THR  %10.4g%10.1f%10.1f\n", -cut, 32.0, 32.0);
-       // 34.0 = Xi zero
-       // 35.0 = AntiXi zero
-       fprintf(fgFile,"PART-THR  %10.4g%10.1f%10.1f\n", -cut, 34.0, 35.0);
-       // 47.0 = D zero
-       // 48.0 = AntiD zero
-       fprintf(fgFile,"PART-THR  %10.4g%10.1f%10.1f\n", -cut, 47.0, 48.0);
-       // 53.0 = Xi_c zero
-       fprintf(fgFile,"PART-THR  %10.4g%10.1f%10.1f\n", -cut, 53.0, 53.0);
-       // 55.0 = Xi'_c zero
-       // 56.0 = Omega_c zero
-       fprintf(fgFile,"PART-THR  %10.4g%10.1f%10.1f\n", -cut, 55.0, 56.0);
-       // 59.0 = AntiXi_c zero
-       fprintf(fgFile,"PART-THR  %10.4g%10.1f%10.1f\n", -cut, 59.0, 59.0);
-       // 61.0 = AntiXi'_c zero
-       // 62.0 = AntiOmega_c zero
-       fprintf(fgFile,"PART-THR  %10.4g%10.1f%10.1f\n", -cut, 61.0, 62.0);
+  // Cut on neutral hadrons
+  fprintf(fgFile,"*\n*Cut for neutral hadrons. CUTNEU = %13.4g\n", fCutValue[kCUTNEU]);
+  
+  // Cut on neutral hadrons
+  fprintf(fgFile,"*\n*Cut for neutral hadrons. CUTNEU = %13.4g\n", fCutValue[kCUTNEU]);
+  
+  // Energy group structure of the 72-neutron ENEA data set:
+  const Float_t neutronGroupUpLimit[72] = {
+    1.9600E-02, 1.7500E-02, 1.4918E-02, 1.3499E-02, 1.2214E-02, 1.1052E-02, 1.0000E-02, 9.0484E-03,
+    8.1873E-03, 7.4082E-03, 6.7032E-03, 6.0653E-03, 5.4881E-03, 4.9659E-03, 4.4933E-03, 4.0657E-03,
+    3.6788E-03, 3.3287E-03, 3.0119E-03, 2.7253E-03, 2.4660E-03, 2.2313E-03, 2.0190E-03, 1.8268E-03,
+    1.6530E-03, 1.4957E-03, 1.3534E-03, 1.2246E-03, 1.1080E-03, 1.0026E-03, 9.0718E-04, 8.2085E-04,
+    7.4274E-04, 6.0810E-04, 4.9787E-04, 4.0762E-04, 3.3373E-04, 2.7324E-04, 2.2371E-04, 1.8316E-04,
+    1.4996E-04, 1.2277E-04, 8.6517E-05, 5.2475E-05, 3.1828E-05, 2.1852E-05, 1.5034E-05, 1.0332E-05,
+    7.1018E-06, 4.8809E-06, 3.3546E-06, 2.3054E-06, 1.5846E-06, 1.0446E-06, 6.8871E-07, 4.5400E-07,
+    2.7537E-07, 1.6702E-07, 1.0130E-07, 6.1442E-08, 3.7267E-08, 2.2603E-08, 1.5535E-08, 1.0677E-08,
+    7.3375E-09, 5.0435E-09, 3.4662E-09, 2.3824E-09, 1.6374E-09, 1.1254E-09, 6.8257E-10, 4.1400E-10
+  };
+
+  Float_t cut = fCutValue[kCUTNEU];
+  //
+  // 8.0 = Neutron
+  // 9.0 = Antineutron
+  // Find the FLUKA neutron group corresponding to the cut
+  //
+  Float_t neutronCut = cut;
+  Int_t groupCut = 1; // if cut is > 19.6 MeV not low energy neutron transport is done
+  if (neutronCut < 0.0196) {
+    neutronCut = 0.0196;
+    // Search the group cutoff for the energy cut
+    Int_t i;
+    for( i=0; i<72; i++ ) {
+      if (cut > neutronGroupUpLimit[i]) break;
+    }
+    groupCut = i+1;
+  } 
+  
+  if (fMedium == -1) {
+        Float_t cut = fCutValue[kCUTNEU];
+        // 8.0 = Neutron
+        // 9.0 = Antineutron
+        fprintf(fgFile,"PART-THR  %10.4g%10.1f%10.1f\n", -neutronCut,  8.0,  9.0);
+        fprintf(fgFile,"LOW-BIAS  %10.4g%10.4g%10.1f%10.1f%10.1f%10.1f\n",
+              Float_t(groupCut), 73.0, 0.95, 2., Float_t(fgGeom->NofVolumes()), 1.);
+        //
+       //
+       // 12.0 = Kaon zero long
+       fprintf(fgFile,"PART-THR  %10.4g%10.1f%10.1f\n", -cut, 12.0, 12.0);
+       // 17.0 = Lambda, 18.0 = Antilambda
+       // 19.0 = Kaon zero short
+       fprintf(fgFile,"PART-THR  %10.4g%10.1f%10.1f\n", -cut, 17.0, 19.0);
+       // 22.0 = Sigma zero, Pion zero, Kaon zero
+       // 25.0 = Antikaon zero
+       fprintf(fgFile,"PART-THR  %10.4g%10.1f%10.1f\n", -cut, 22.0, 25.0);
+       // 32.0 = Antisigma zero
+       fprintf(fgFile,"PART-THR  %10.4g%10.1f%10.1f\n", -cut, 32.0, 32.0);
+       // 34.0 = Xi zero
+       // 35.0 = AntiXi zero
+       fprintf(fgFile,"PART-THR  %10.4g%10.1f%10.1f\n", -cut, 34.0, 35.0);
+       // 47.0 = D zero
+       // 48.0 = AntiD zero
+       fprintf(fgFile,"PART-THR  %10.4g%10.1f%10.1f\n", -cut, 47.0, 48.0);
+       // 53.0 = Xi_c zero
+       fprintf(fgFile,"PART-THR  %10.4g%10.1f%10.1f\n", -cut, 53.0, 53.0);
+       // 55.0 = Xi'_c zero
+       // 56.0 = Omega_c zero
+       fprintf(fgFile,"PART-THR  %10.4g%10.1f%10.1f\n", -cut, 55.0, 56.0);
+       // 59.0 = AntiXi_c zero
+       fprintf(fgFile,"PART-THR  %10.4g%10.1f%10.1f\n", -cut, 59.0, 59.0);
+       // 61.0 = AntiXi'_c zero
+       // 62.0 = AntiOmega_c zero
+       fprintf(fgFile,"PART-THR  %10.4g%10.1f%10.1f\n", -cut, 61.0, 62.0);
     } else {
-       printf("Cuts on neutral hadrons material by material not yet implemented !\n");
+        Int_t nreg, *reglist;
+        Float_t ireg;
+        reglist = fgGeom->GetMaterialList(fMedium, nreg);
+        // Loop over regions of a given material
+        for (Int_t k = 0; k < nreg; k++) {
+         ireg = reglist[k];
+         fprintf(fgFile,"LOW-BIAS  %10.4g%10.4g%10.1f%10.1f%10.1f%10.1f\n",
+                Float_t(groupCut), 73.0, 0.95, ireg, ireg, 1.);
+       }
+
+       Warning("ProcessCUTNEU",
+              "Material #%4d %s: Cut on neutral hadrons (Ekin > %9.3e) material by material only implemented for low-energy neutrons !\n",
+              fMedium, fCMaterial->GetName(), cut);
     }
 }
 
@@ -593,39 +755,41 @@ void TFlukaConfigOption::ProcessCUTHAD()
 {
     // Cut on charged hadrons
     fprintf(fgFile,"*\n*Cut for charge hadrons. CUTHAD = %13.4g\n", fCutValue[kCUTHAD]);
+    Float_t cut = fCutValue[kCUTHAD];
     if (fMedium == -1) {
-       Float_t cut = fCutValue[kCUTHAD];
-       // 1.0 = Proton
-       // 2.0 = Antiproton
-       fprintf(fgFile,"PART-THR  %10.4g%10.1f%10.1f\n", -cut,  1.0,  2.0);
-       // 13.0 = Positive Pion, Negative Pion, Positive Kaon
-       // 16.0 = Negative Kaon
-       fprintf(fgFile,"PART-THR  %10.4g%10.1f%10.1f\n", -cut, 13.0, 16.0);
-       // 20.0 = Negative Sigma
-       // 21.0 = Positive Sigma
-       fprintf(fgFile,"PART-THR  %10.4g%10.1f%10.1f\n", -cut, 20.0, 21.0);
-       // 31.0 = Antisigma minus
-       // 33.0 = Antisigma plus
-       fprintf(fgFile,"PART-THR  %10.4g%10.1f%10.1f\n", -cut, 31.0, 31.0);
-       fprintf(fgFile,"PART-THR  %10.4g%10.1f%10.1f\n", -cut, 33.0, 33.0);
-       // 36.0 = Negative Xi, Positive Xi, Omega minus
-       // 39.0 = Antiomega
-       fprintf(fgFile,"PART-THR  %10.4g%10.1f%10.1f\n", -cut, 36.0, 39.0);
-       // 45.0 = D plus
-       // 46.0 = D minus
-       fprintf(fgFile,"PART-THR  %10.4g%10.1f%10.1f\n", -cut, 45.0, 46.0);
-       // 49.0 = D_s plus, D_s minus, Lambda_c plus
-       // 52.0 = Xi_c plus
-       fprintf(fgFile,"PART-THR  %10.4g%10.1f%10.1f\n", -cut, 49.0, 52.0);
-       // 54.0 = Xi'_c plus
-       // 60.0 = AntiXi'_c minus
-       fprintf(fgFile,"PART-THR  %10.4g%10.1f%10.1f\n", -cut, 54.0, 54.0);
-       fprintf(fgFile,"PART-THR  %10.4g%10.1f%10.1f\n", -cut, 60.0, 60.0);
-       // 57.0 = Antilambda_c minus
-       // 58.0 = AntiXi_c minus
-       fprintf(fgFile,"PART-THR  %10.4g%10.1f%10.1f\n", -cut, 57.0, 58.0);
+       // 1.0 = Proton
+       // 2.0 = Antiproton
+       fprintf(fgFile,"PART-THR  %10.4g%10.1f%10.1f\n", -cut,  1.0,  2.0);
+       // 13.0 = Positive Pion, Negative Pion, Positive Kaon
+       // 16.0 = Negative Kaon
+       fprintf(fgFile,"PART-THR  %10.4g%10.1f%10.1f\n", -cut, 13.0, 16.0);
+       // 20.0 = Negative Sigma
+       // 21.0 = Positive Sigma
+       fprintf(fgFile,"PART-THR  %10.4g%10.1f%10.1f\n", -cut, 20.0, 21.0);
+       // 31.0 = Antisigma minus
+       // 33.0 = Antisigma plus
+       fprintf(fgFile,"PART-THR  %10.4g%10.1f%10.1f\n", -cut, 31.0, 31.0);
+       fprintf(fgFile,"PART-THR  %10.4g%10.1f%10.1f\n", -cut, 33.0, 33.0);
+       // 36.0 = Negative Xi, Positive Xi, Omega minus
+       // 39.0 = Antiomega
+       fprintf(fgFile,"PART-THR  %10.4g%10.1f%10.1f\n", -cut, 36.0, 39.0);
+       // 45.0 = D plus
+       // 46.0 = D minus
+       fprintf(fgFile,"PART-THR  %10.4g%10.1f%10.1f\n", -cut, 45.0, 46.0);
+       // 49.0 = D_s plus, D_s minus, Lambda_c plus
+       // 52.0 = Xi_c plus
+       fprintf(fgFile,"PART-THR  %10.4g%10.1f%10.1f\n", -cut, 49.0, 52.0);
+       // 54.0 = Xi'_c plus
+       // 60.0 = AntiXi'_c minus
+       fprintf(fgFile,"PART-THR  %10.4g%10.1f%10.1f\n", -cut, 54.0, 54.0);
+       fprintf(fgFile,"PART-THR  %10.4g%10.1f%10.1f\n", -cut, 60.0, 60.0);
+       // 57.0 = Antilambda_c minus
+       // 58.0 = AntiXi_c minus
+       fprintf(fgFile,"PART-THR  %10.4g%10.1f%10.1f\n", -cut, 57.0, 58.0);
     } else {
-       printf("Cuts on charged hadrons material by material not yet implemented !\n");
+      Warning("ProcessCUTHAD", 
+              "Material #%4d %s: Cut on charged hadrons (Ekin > %9.3e) material by material not yet implemented !\n",
+             fMedium, fCMaterial->GetName(), cut);
     }
 }
 
@@ -635,9 +799,10 @@ void TFlukaConfigOption::ProcessCUTMUO()
     fprintf(fgFile,"*\n*Cut for muons. CUTMUO = %13.4g\n", fCutValue[kCUTMUO]);
     Float_t cut = fCutValue[kCUTMUO];
     if (fMedium == -1) {
-       fprintf(fgFile,"PART-THR  %10.4g%10.1f%10.1f\n",-cut, 10.0, 11.0);
+       fprintf(fgFile,"PART-THR  %10.4g%10.1f%10.1f\n",-cut, 10.0, 11.0);
     } else {
-       printf("Cuts on muons material by material not yet implemented !\n");
+       Warning("ProcessCUTMUO", "Material #%4d %s: Cut on muons (Ekin > %9.3e) material by material not yet implemented !\n",
+              fMedium, fCMaterial->GetName(), cut);
     }
     
     
@@ -650,3 +815,18 @@ void TFlukaConfigOption::ProcessTOFMAX()
     fprintf(fgFile,"*\n*Cut on time of flight. TOFMAX = %13.4g\n", fCutValue[kTOFMAX]);
     fprintf(fgFile,"TIME-CUT  %10.4g%10.1f%10.1f%10.1f%10.1f\n",cut*1.e9,0.,0.,-6.0,64.0);
 }
+
+void  TFlukaConfigOption::ProcessSensitiveMedium()
+{
+    //
+    // Special options for sensitive media
+    //
+
+    fprintf(fgFile,"*\n*Options for sensitive medium \n");
+    //
+    // EMFFIX
+    fprintf(fgFile,"EMFFIX    %10.4g%10.4g%10.1f%10.1f%10.1f%10.1f\n", fCMatMin, 0.05, 0., 0., 0., 0.);
+    //
+    // FLUKAFIX
+    fprintf(fgFile,"FLUKAFIX  %10.3g                    %10.3f\n", 0.05, fCMatMin);
+}