Writing of FLUKA input cards for physics configuration is now the responsibility of
authormorsch <morsch@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 16 Dec 2004 09:42:39 +0000 (09:42 +0000)
committermorsch <morsch@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 16 Dec 2004 09:42:39 +0000 (09:42 +0000)
TFlukaConfigOptions. One TFlukaConfigOption per material with special process flags and cuts
is created. One object with material id = -1 is created to hold the global settings.

TFluka/TFluka.cxx
TFluka/TFluka.h
TFluka/TFlukaConfigOption.cxx
TFluka/TFlukaConfigOption.h

index 0f48827c017ba86e17aa1be415bbad54929bcfc5..261289ca38da8a910313f6ef817d8032f9ecd19f 100644 (file)
@@ -103,8 +103,7 @@ TFluka::TFluka()
   :TVirtualMC(),
    fVerbosityLevel(0),
    fInputFileName(""),
-   fProcesses(0), 
-   fCuts(0),
+   fUserConfig(0), 
    fUserScore(0)
 { 
   //
@@ -132,8 +131,7 @@ TFluka::TFluka(const char *title, Int_t verbosity, Bool_t isRootGeometrySupporte
    fTrackIsEntering(0),
    fTrackIsExiting(0),
    fTrackIsNew(0),
-   fProcesses(new TObjArray(100)),
-   fCuts(new TObjArray(100)), 
+   fUserConfig(new TObjArray(100)),
    fUserScore(new TObjArray(100)) 
 {
   // create geometry interface
@@ -148,7 +146,7 @@ TFluka::TFluka(const char *title, Int_t verbosity, Bool_t isRootGeometrySupporte
    fFieldFlag = 1;
    fGeneratePemf = kFALSE;
    fMCGeo = new TGeoMCGeometry("MCGeo", "TGeo Implementation of VirtualMCGeometry", kTRUE);
-   fGeom = new TFlukaMCGeometry("geom", "FLUKA VMC Geometry");
+   fGeom  = new TFlukaMCGeometry("geom", "FLUKA VMC Geometry");
    if (verbosity > 2) fGeom->SetDebugMode(kTRUE);
    fMaterials = 0;
    fStopped   = 0;
@@ -166,14 +164,9 @@ TFluka::~TFluka() {
     delete fGeom;
     delete fMCGeo;
     
-    if (fCuts) {
-       fCuts->Delete();
-       delete fCuts;
-    }
-
-    if (fProcesses) {
-       fProcesses->Delete();
-       delete fProcesses;
+    if (fUserConfig) {
+       fUserConfig->Delete();
+       delete fUserConfig;
     }
 
 
@@ -378,6 +371,7 @@ void TFluka::Material(Int_t& kmat, const char* name, Double_t a,
                      Double_t z, Double_t dens, Double_t radl, Double_t absl,
                      Double_t* /*buf*/, Int_t /*nwbuf*/) {
 //
+// Define a material
   TGeoMaterial *mat;
   kmat = gGeoManager->GetListOfMaterials()->GetSize();
   if ((z-Int_t(z)) > 1E-3) {
@@ -394,6 +388,8 @@ void TFluka::Material(Int_t& kmat, const char* name, Double_t a,
 //______________________________________________________________________________ 
 void TFluka::Mixture(Int_t& kmat, const char *name, Float_t *a, 
                     Float_t *z, Double_t dens, Int_t nlmat, Float_t *wmat) {
+//
+// Define a material mixture
 //
   Double_t* da = fGeom->CreateDoubleArray(a, TMath::Abs(nlmat));  
   Double_t* dz = fGeom->CreateDoubleArray(z, TMath::Abs(nlmat));  
@@ -819,37 +815,31 @@ void TFluka::StopTrack()
 void TFluka::SetProcess(const char* flagName, Int_t flagValue, Int_t imed)
 {
 //  Set process user flag for material imat
-//
-    TFlukaConfigOption* proc = new TFlukaConfigOption(flagName, flagValue, imed);
-    fProcesses->Add(proc);
-}
-
-//______________________________________________________________________________ 
-Bool_t TFluka::SetProcess(const char* flagName, Int_t flagValue)
-{
-//  Set process user flag 
 //
 //    
 //  Update if already in the list
 //
-
-    TIter next(fProcesses);
+    TIter next(fUserConfig);
     TFlukaConfigOption* proc;
     while((proc = (TFlukaConfigOption*)next()))
     { 
-       if (strcmp(proc->GetName(), flagName) == 0) {
-           proc->SetFlag(flagValue);
-           proc->SetMedium(-1);
-           return kTRUE;
-         }
+       if (proc->Medium() == imed) {
+           proc->SetProcess(flagName, flagValue);
+           return;
+       }
     }
+    proc = new TFlukaConfigOption(imed);
+    proc->SetProcess(flagName, flagValue);
+    fUserConfig->Add(proc);
+}
+
+//______________________________________________________________________________ 
+Bool_t TFluka::SetProcess(const char* flagName, Int_t flagValue)
+{
+//  Set process user flag 
 //
-// If not create a new process
 //    
-
-    proc = new TFlukaConfigOption(flagName, flagValue);
-    fProcesses->Add(proc);
-    
+    SetProcess(flagName, flagValue, -1);
     return kTRUE;  
 }
 
@@ -858,8 +848,19 @@ void TFluka::SetCut(const char* cutName, Double_t cutValue, Int_t imed)
 {
 // Set user cut value for material imed
 //
-    TFlukaConfigOption* cut = new TFlukaConfigOption(cutName, cutValue, imed);
-    fCuts->Add(cut);
+    TIter next(fUserConfig);
+    TFlukaConfigOption* proc;
+    while((proc = (TFlukaConfigOption*)next()))
+    { 
+       if (proc->Medium() == imed) {
+           proc->SetCut(cutName, cutValue);
+           return;
+       }
+    }
+
+    proc = new TFlukaConfigOption(imed);
+    proc->SetCut(cutName, cutValue);
+    fUserConfig->Add(proc);
 }
 
 //______________________________________________________________________________ 
@@ -868,26 +869,8 @@ Bool_t TFluka::SetCut(const char* cutName, Double_t cutValue)
 // Set user cut value 
 //
 //    
-//  Update if already in the list
-//
-
-    TIter next(fCuts);
-    TFlukaConfigOption* cut;
-    while((cut = (TFlukaConfigOption*)next()))
-    { 
-       if (strcmp(cut->GetName(), cutName) == 0) {
-           cut->SetCut(cutValue);
-           return kTRUE;
-         }
-    }
-//
-// If not create a new process
-//    
-
-    cut = new TFlukaConfigOption(cutName, cutValue);
-    fCuts->Add(cut);
-    
-    return kTRUE;  
+    SetCut(cutName, cutValue, -1);
+    return kTRUE;
 }
 
 void TFluka::SetUserScoring(const char* option, Int_t npar, Float_t what[12])
@@ -913,1117 +896,92 @@ void TFluka::InitPhysics()
 //
 // Physics initialisation with preparation of FLUKA input cards
 //
-  printf("=>InitPhysics\n");
-  Int_t j, k;
-  Double_t theCut;
-
-  FILE *pFlukaVmcCoreInp, *pFlukaVmcFlukaMat, *pFlukaVmcInp;
-
-  Double_t zero  = 0.0;
-  Double_t one   = 1.0;
-  Double_t two   = 2.0;
-  Double_t three = 3.0;
-
-  Float_t fLastMaterial = fGeom->GetLastMaterialIndex();
-  if (fVerbosityLevel >= 3) printf("   last FLUKA material is %g\n", fLastMaterial);
-
-  // Prepare  Cerenkov
-  TObjArray *matList = GetFlukaMaterials();
-  Int_t nmaterial =  matList->GetEntriesFast();
-  fMaterials = new Int_t[nmaterial+3];
-             
-// construct file names
+    printf("=>InitPhysics\n");
 
-  TString sFlukaVmcCoreInp = getenv("ALICE_ROOT");
-  sFlukaVmcCoreInp +="/TFluka/input/";
-  TString sFlukaVmcTmp = "flukaMat.inp";
-  TString sFlukaVmcInp = GetInputFileName();
-  sFlukaVmcCoreInp += GetCoreInputFileName();
-
-// open files 
-
-  if ((pFlukaVmcCoreInp = fopen(sFlukaVmcCoreInp.Data(),"r")) == NULL) {
-      printf("\nCannot open file %s\n",sFlukaVmcCoreInp.Data());
-      exit(1);
-  }
-  if ((pFlukaVmcFlukaMat = fopen(sFlukaVmcTmp.Data(),"r")) == NULL) {
-      printf("\nCannot open file %s\n",sFlukaVmcTmp.Data());
-      exit(1);
-  }
-  if ((pFlukaVmcInp = fopen(sFlukaVmcInp.Data(),"w")) == NULL) {
-      printf("\nCannot open file %s\n",sFlukaVmcInp.Data());
-      exit(1);
-  }
+// Construct file names
+    FILE *pFlukaVmcCoreInp, *pFlukaVmcFlukaMat, *pFlukaVmcInp;
+    TString sFlukaVmcCoreInp = getenv("ALICE_ROOT");
+    sFlukaVmcCoreInp +="/TFluka/input/";
+    TString sFlukaVmcTmp = "flukaMat.inp";
+    TString sFlukaVmcInp = GetInputFileName();
+    sFlukaVmcCoreInp += GetCoreInputFileName();
+    
+// Open files 
+    if ((pFlukaVmcCoreInp = fopen(sFlukaVmcCoreInp.Data(),"r")) == NULL) {
+       printf("\nCannot open file %s\n",sFlukaVmcCoreInp.Data());
+       exit(1);
+    }
+    if ((pFlukaVmcFlukaMat = fopen(sFlukaVmcTmp.Data(),"r")) == NULL) {
+       printf("\nCannot open file %s\n",sFlukaVmcTmp.Data());
+       exit(1);
+    }
+    if ((pFlukaVmcInp = fopen(sFlukaVmcInp.Data(),"w")) == NULL) {
+       printf("\nCannot open file %s\n",sFlukaVmcInp.Data());
+       exit(1);
+    }
 
-// copy core input file 
-  Char_t sLine[255];
-  Float_t fEventsPerRun;
-  
-  while ((fgets(sLine,255,pFlukaVmcCoreInp)) != NULL) {
-      if (strncmp(sLine,"GEOEND",6) != 0)
-         fprintf(pFlukaVmcInp,"%s",sLine); // copy until GEOEND card
-      else {
-         fprintf(pFlukaVmcInp,"GEOEND\n");   // add GEOEND card
-         goto flukamat;
-      }
-  } // end of while until GEOEND card
-  
+// Copy core input file 
+    Char_t sLine[255];
+    Float_t fEventsPerRun;
+    
+    while ((fgets(sLine,255,pFlukaVmcCoreInp)) != NULL) {
+       if (strncmp(sLine,"GEOEND",6) != 0)
+           fprintf(pFlukaVmcInp,"%s",sLine); // copy until GEOEND card
+       else {
+           fprintf(pFlukaVmcInp,"GEOEND\n");   // add GEOEND card
+           goto flukamat;
+       }
+    } // end of while until GEOEND card
+    
 
  flukamat:
-  while ((fgets(sLine,255,pFlukaVmcFlukaMat)) != NULL) { // copy flukaMat.inp file
-      fprintf(pFlukaVmcInp,"%s\n",sLine);
-  }
-  
-  while ((fgets(sLine,255,pFlukaVmcCoreInp)) != NULL) { 
-      if (strncmp(sLine,"START",5) != 0)
-         fprintf(pFlukaVmcInp,"%s\n",sLine);
-      else {
-         sscanf(sLine+10,"%10f",&fEventsPerRun);
-      goto fin;
-      }
-  } //end of while until START card
-  
-fin:
-// in G3 the process control values meaning can be different for
-// different processes, but for most of them is:
-//   0  process is not activated
-//   1  process is activated WITH generation of secondaries
-//   2  process is activated WITHOUT generation of secondaries
-// if process does not generate secondaries => 1 same as 2
-//
-// Exceptions:
-//   MULS:  also 3
-//   LOSS:  also 3, 4
-//   RAYL:  only 0,1
-//   HADR:  may be > 2
-//
-// Loop over number of SetProcess calls 
-  fprintf(pFlukaVmcInp,"*----------------------------------------------------------------------------- \n");
-  fprintf(pFlukaVmcInp,"*----- The following data are generated from SetProcess and SetCut calls ----- \n");
-  fprintf(pFlukaVmcInp,"*----------------------------------------------------------------------------- \n");
-
-// Outer loop over processes
-  TIter next(fProcesses);
-  TFlukaConfigOption *proc;
-// Inner loop over processes
-  TIter nextp(fProcesses);
-  TFlukaConfigOption *procp;
-// Loop over cuts
-  TIter nextc(fCuts);
-  TFlukaConfigOption *cut = 0x0;
-
-  while((proc = (TFlukaConfigOption*)next())) {
-      Float_t matMin = three;
-      Float_t matMax = fLastMaterial;
-      Bool_t  global = kTRUE;
-      if (proc->Medium() != -1) {
-         Int_t mat;
-         if ((mat = proc->Medium()) >= GetFlukaMaterials()->GetEntries()) continue;
-         matMin = Float_t(mat);
-         matMax = matMin;
-         global = kFALSE;
-
-         fprintf(pFlukaVmcInp,"*\n*Material specific process setting for #%8d \n", mat);
-      }
-      
-    // annihilation
-    // G3 default value: 1
-    // G4 processes: G4eplusAnnihilation/G4IeplusAnnihilation
-    // Particles: e+
-    // Physics:   EM
-    // flag = 0 no annihilation
-    // flag = 1 annihilation, decays processed
-    // flag = 2 annihilation, no decay product stored
-    // gMC ->SetProcess("ANNI",1); // EMFCUT   -1.   0.  0. 3. lastmat 0. ANNH-THR
-      if (strncmp(proc->GetName(),"ANNI",4) == 0) {
-         if (proc->Flag() == 1 || proc->Flag() == 2) {
-             fprintf(pFlukaVmcInp,"*\n*Kinetic energy threshold (GeV) for e+ annihilation - resets to default=0.\n");
-             fprintf(pFlukaVmcInp,"*Generated from call: SetProcess('ANNI',1) or SetProcess('ANNI',2)\n");
-             // -one = kinetic energy threshold (GeV) for e+ annihilation (resets to default=0)
-             // zero = not used
-             // zero = not used
-             // matMin = lower bound of the material indices in which the respective thresholds apply
-             // matMax = upper bound of the material indices in which the respective thresholds apply
-             // one = step length in assigning indices
-             // "ANNH-THR"; 
-             fprintf(pFlukaVmcInp,"EMFCUT    %10.1f%10.1f%10.1f%10.1f%10.1f%10.1fANNH-THR\n",-one,zero,zero,matMin,matMax,one);
-         }
-         else if (proc->Flag() == 0) {
-             fprintf(pFlukaVmcInp,"*\n*No annihilation - no FLUKA card generated\n");
-             fprintf(pFlukaVmcInp,"*Generated from call: SetProcess('ANNI',0)\n");
-         }
-         else  {
-             fprintf(pFlukaVmcInp,"*\n*Illegal flag value in SetProcess('ANNI',?) call.\n");
-             fprintf(pFlukaVmcInp,"*No FLUKA card generated\n");
-         }
-      }
+    while ((fgets(sLine,255,pFlukaVmcFlukaMat)) != NULL) { // copy flukaMat.inp file
+       fprintf(pFlukaVmcInp,"%s\n",sLine);
+    }
     
-    // bremsstrahlung and pair production are both activated
-    // G3 default value: 1
-    // G4 processes: G4eBremsstrahlung/G4IeBremsstrahlung,
-    //               G4MuBremsstrahlung/G4IMuBremsstrahlung,
-    //               G4LowEnergyBremstrahlung
-    // Particles: e-/e+; mu+/mu-
-    // Physics:   EM
-    // flag = 0 no bremsstrahlung
-    // flag = 1 bremsstrahlung, photon processed
-    // flag = 2 bremsstrahlung, no photon stored
-    // gMC ->SetProcess("BREM",1); // PAIRBREM  2.   0.  0. 3. lastmat
-                                 // EMFCUT   -1.   0.  0. 3. lastmat 0. ELPO-THR
-    // G3 default value: 1
-    // G4 processes: G4GammaConversion,
-    //               G4MuPairProduction/G4IMuPairProduction
-    //               G4LowEnergyGammaConversion
-    // Particles: gamma, mu
-    // Physics:   EM
-    // flag = 0 no delta rays
-    // flag = 1 delta rays, secondaries processed
-    // flag = 2 delta rays, no secondaries stored
-
-    else if ((strncmp(proc->GetName(),"PAIR",4) == 0) && (proc->Flag() == 1 || proc->Flag() == 2)) {
-
-       nextp.Reset();
-       
-       while ((procp = (TFlukaConfigOption*)nextp())) {
-           if ((strncmp(procp->GetName(),"BREM",4) == 0) && 
-               (proc->Flag() == 1) &&
-               (procp->Medium() == proc->Medium())) {
-               fprintf(pFlukaVmcInp,"*\n*Bremsstrahlung and pair production by muons and charged hadrons both activated\n");
-               fprintf(pFlukaVmcInp,"*Generated from call: SetProcess('BREM',1) and SetProcess('PAIR',1)\n");
-               fprintf(pFlukaVmcInp,"*Energy threshold set by call SetCut('BCUTM',cut) or set to 0.\n");
-               fprintf(pFlukaVmcInp,"*Energy threshold set by call SetCut('PPCUTM',cut) or set to 0.\n");
-               // three = bremsstrahlung and pair production by muons and charged hadrons both are activated
-               fprintf(pFlukaVmcInp,"PAIRBREM  %10.1f",three);
-               // direct pair production by muons
-               // G4 particles: "e-", "e+"
-               // G3 default value: 0.01 GeV
-               //gMC ->SetCut("PPCUTM",cut); // total energy cut for direct pair prod. by muons
-               theCut = 0.0;
-               nextc.Reset();
-               while ((cut = (TFlukaConfigOption*)nextc())) {
-                   if (strncmp(cut->GetName(), "PPCUTM", 6) == 0 &&
-                       (cut->Medium() == proc->Medium())) theCut = cut->Cut();
-               }
-               fprintf(pFlukaVmcInp,"%10.4g",theCut);
-               // theCut; = e+, e- kinetic energy threshold (in GeV) for explicit pair production.
-               // muon and hadron bremsstrahlung
-               // G4 particles: "gamma"
-               // G3 default value: CUTGAM=0.001 GeV
-               //gMC ->SetCut("BCUTM",cut);  // cut for muon and hadron bremsstrahlung
-               theCut = 0.0;
-               nextc.Reset();
-               while ((cut = (TFlukaConfigOption*)nextc())) {
-                   if (strncmp(cut->GetName(), "BCUTM", 5) == 0 &&
-                       (cut->Medium() == proc->Medium())) theCut = cut->Cut();
-               }
-               fprintf(pFlukaVmcInp,"%10.4g%10.1f%10.1f\n",theCut,matMin,matMax);
-               // theCut = photon energy threshold (GeV) for explicit bremsstrahlung production
-               // matMin = lower bound of the material indices in which the respective thresholds apply
-               // matMax = upper bound of the material indices in which the respective thresholds apply
-               
-               // for e+ and e-
-               fprintf(pFlukaVmcInp,"*\n*Kinetic energy threshold (GeV) for e+/e- bremsstrahlung - resets to default=0.\n");
-               fprintf(pFlukaVmcInp,"*Generated from call: SetProcess('BREM',1);\n");
-               theCut = -1.0;
-               nextc.Reset();
-               while ((cut = (TFlukaConfigOption*)nextc())) {
-                   if (strncmp(cut->GetName(), "BCUTE", 5) == 0 &&
-                       (cut->Medium() == proc->Medium())) theCut = cut->Cut();
-               }
-               //theCut = kinetic energy threshold (GeV) for e+/e- bremsstrahlung (resets to default=0)
-               // zero = not used
-               // zero = not used
-               // matMin = lower bound of the material indices in which the respective thresholds apply
-               // matMax = upper bound of the material indices in which the respective thresholds apply
-               // one = step length in assigning indices
-               // "ELPO-THR"; 
-               fprintf(pFlukaVmcInp,"EMFCUT    %10.4g%10.1f%10.1f%10.1f%10.1f%10.1fELPO-THR\n",theCut,zero,zero,matMin,matMax,one);
-               
-                // for gamma -> e+ and e-
-               fprintf(pFlukaVmcInp,"*\n*Pair production by photons is activated\n");
-               fprintf(pFlukaVmcInp,"*Generated from call: SetProcess('PAIR',1);\n");
-               theCut = -1.0;
-               nextc.Reset();
-               while ((cut = (TFlukaConfigOption*)nextc())) {
-                   if (strncmp(cut->GetName(), "CUTELE", 6) == 0 &&
-                       (cut->Medium() == proc->Medium())) theCut = cut->Cut();
-               }
-               // theCut = energy threshold (GeV) for gamma pair production (< 0.0 : resets to default, = 0.0 : ignored)
-               // matMin = lower bound of the material indices in which the respective thresholds apply
-               // matMax =  upper bound of the material indices in which the respective thresholds apply
-               // one = step length in assigning indices
-               fprintf(pFlukaVmcInp,"EMFCUT    %10.1f%10.1f%10.4g%10.1f%10.1f%10.1fPHOT-THR\n",zero,zero,theCut,matMin,matMax,one);
-               goto BOTH;
-           } // end of if for BREM
-       } // end of loop for BREM
-       
-       // only pair production by muons and charged hadrons is activated
-       fprintf(pFlukaVmcInp,"*\n*Pair production by muons and charged hadrons is activated\n");
-       fprintf(pFlukaVmcInp,"*Generated from call: SetProcess('PAIR',1) or SetProcess('PAIR',2)\n");
-       fprintf(pFlukaVmcInp,"*Energy threshold set by call SetCut('PPCUTM',cut) or set to 0.\n");
-       // direct pair production by muons
-       // G4 particles: "e-", "e+"
-       // G3 default value: 0.01 GeV
-       //gMC ->SetCut("PPCUTM",cut); // total energy cut for direct pair prod. by muons
-       // one = pair production by muons and charged hadrons is activated
-       // zero = e+, e- kinetic energy threshold (in GeV) for explicit pair production.
-       // zero = no explicit bremsstrahlung production is simulated
-       // matMin = lower bound of the material indices in which the respective thresholds apply
-       // matMax = upper bound of the material indices in which the respective thresholds apply
-       fprintf(pFlukaVmcInp,"PAIRBREM  %10.1f%10.1f%10.1f%10.1f%10.1f\n",one,zero,zero,matMin,matMax);
-       
-       // for gamma -> e+ and e-
-       fprintf(pFlukaVmcInp,"*\n*Pair production by electrons is activated\n");
-       fprintf(pFlukaVmcInp,"*Generated from call: SetProcess('PAIR',1) or SetProcess('PAIR',2)\n");
-       theCut = -1.0;
-       nextc.Reset();
-       while ((cut = (TFlukaConfigOption*)nextc())) {
-           if (strncmp(cut->GetName(), "CUTELE", 6) == 0 &&
-               (cut->Medium() == proc->Medium())) theCut = cut->Cut();
+    while ((fgets(sLine,255,pFlukaVmcCoreInp)) != NULL) { 
+       if (strncmp(sLine,"START",5) != 0)
+           fprintf(pFlukaVmcInp,"%s\n",sLine);
+       else {
+           sscanf(sLine+10,"%10f",&fEventsPerRun);
+           goto fin;
        }
-       // zero = energy threshold (GeV) for Compton scattering (= 0.0 : ignored)
-       // zero = energy threshold (GeV) for Photoelectric (= 0.0 : ignored)
-       // theCut = energy threshold (GeV) for gamma pair production (< 0.0 : resets to default, = 0.0 : ignored)
-       // matMin = lower bound of the material indices in which the respective thresholds apply
-       // matMax = upper bound of the material indices in which the respective thresholds apply
-       // one = step length in assigning indices
-       fprintf(pFlukaVmcInp,"EMFCUT    %10.1f%10.1f%10.4g%10.1f%10.1f%10.1fPHOT-THR\n",zero,zero,theCut,matMin,matMax,one);
-      
-    BOTH:
-       k = 0;
-    } // end of if for PAIR
-      
-      
-      
-      // bremsstrahlung
-      // G3 default value: 1
-      // G4 processes: G4eBremsstrahlung/G4IeBremsstrahlung,
-      //               G4MuBremsstrahlung/G4IMuBremsstrahlung,
-      //               G4LowEnergyBremstrahlung
-      // Particles: e-/e+; mu+/mu-
-      // Physics:   EM
-      // flag = 0 no bremsstrahlung
-      // flag = 1 bremsstrahlung, photon processed
-      // flag = 2 bremsstrahlung, no photon stored
-      // gMC ->SetProcess("BREM",1); // PAIRBREM  2.   0.  0. 3. lastmat
-      // EMFCUT   -1.   0.  0. 3. lastmat 0. ELPO-THR
-      else if (strncmp(proc->GetName(),"BREM",4) == 0) {
-         nextp.Reset();
-         while((procp = (TFlukaConfigOption*)nextp())) {
-             if ((strncmp(procp->GetName(),"PAIR",4) == 0) && 
-                 procp->Flag() == 1 &&
-                 (procp->Medium() == proc->Medium())) goto NOBREM;
-         }
-          if (proc->Flag() == 1) { 
-             fprintf(pFlukaVmcInp,"*\n*Bremsstrahlung by muons and charged hadrons is activated\n");
-             fprintf(pFlukaVmcInp,"*Generated from call: SetProcess('BREM',1) or SetProcess('BREM',2)\n");
-             fprintf(pFlukaVmcInp,"*Energy threshold set by call SetCut('BCUTM',cut) or set to 0.\n");
-             // two = bremsstrahlung by muons and charged hadrons is activated
-             // zero = no meaning
-             // muon and hadron bremsstrahlung
-             // G4 particles: "gamma"
-             // G3 default value: CUTGAM=0.001 GeV
-             //gMC ->SetCut("BCUTM",cut);  // cut for muon and hadron bremsstrahlung
-             theCut = 0.0;
-             nextc.Reset();
-             while ((cut = (TFlukaConfigOption*)nextc())) {
-                 if (strncmp(cut->GetName(), "BCUTM", 5) == 0 &&
-                     (cut->Medium() == proc->Medium())) theCut = cut->Cut();
-             }
-             // theCut = photon energy threshold (GeV) for explicit bremsstrahlung production
-             // matMin = lower bound of the material indices in which the respective thresholds apply
-             // matMax = upper bound of the material indices in which the respective thresholds apply
-             fprintf(pFlukaVmcInp,"PAIRBREM  %10.1f%10.1f%10.4g%10.1f%10.1f\n",two,zero,theCut,matMin,matMax);
-             
-             // for e+ and e-
-             fprintf(pFlukaVmcInp,"*\n*Kinetic energy threshold (GeV) for e+/e- bremsstrahlung - resets to default=0.\n");
-             fprintf(pFlukaVmcInp,"*Generated from call: SetProcess('BREM',1);");
-             // - one = kinetic energy threshold (GeV) for e+/e- bremsstrahlung (resets to default=0)
-             // zero = not used
-             // zero = not used
-             // matMin = lower bound of the material indices in which the respective thresholds apply
-             // matMax = upper bound of the material indices in which the respective thresholds apply
-             // one = step length in assigning indices
-             //"ELPO-THR"; 
-             theCut = -1.0;
-             nextc.Reset();
-             while ((cut = (TFlukaConfigOption*)nextc())) {
-                 if (strncmp(cut->GetName(), "CUTGAM", 6) == 0 &&
-                     (cut->Medium() == proc->Medium())) theCut = cut->Cut();
-             }
-             fprintf(pFlukaVmcInp,"EMFCUT    %10.4g%10.1f%10.1f%10.1f%10.1f%10.1fELPO-THR\n", theCut,zero,zero,matMin,matMax,one);
-         }
-         else if (proc->Flag() == 0) {
-             fprintf(pFlukaVmcInp,"*\n*No bremsstrahlung - no FLUKA card generated\n");
-             fprintf(pFlukaVmcInp,"*Generated from call: SetProcess('BREM',0)\n");
-         }
-         else  {
-             fprintf(pFlukaVmcInp,"*\n*Illegal flag value in SetProcess('BREM',?) call.\n");
-             fprintf(pFlukaVmcInp,"*No FLUKA card generated\n");
-         }
-      NOBREM:
-         j = 0;
-      } // end of else if (strncmp(proc->GetName(),"BREM",4) == 0)
-      
-      // Cerenkov photon generation
-      // G3 default value: 0
-      // G4 process: G4Cerenkov
-      // 
-      // Particles: charged
-      // Physics:   Optical
-      // flag = 0 no Cerenkov photon generation
-      // flag = 1 Cerenkov photon generation
-      // flag = 2 Cerenkov photon generation with primary stopped at each step
-      //xx gMC ->SetProcess("CKOV",1); // ??? Cerenkov photon generation
-      
-      else if (strncmp(proc->GetName(),"CKOV",4) == 0) {
-         if ((proc->Flag() == 1 || proc->Flag() == 2) && global) {
-             // Write comments
-             fprintf(pFlukaVmcInp, "* \n"); 
-             fprintf(pFlukaVmcInp, "*Cerenkov photon generation\n"); 
-             fprintf(pFlukaVmcInp, "*Generated from call: SetProcess('CKOV',1) or SetProcess('CKOV',2)\n"); 
-             // Loop over media 
-             for (Int_t im = 0; im < nmaterial; im++)
-             {
-                 TGeoMaterial* material = dynamic_cast<TGeoMaterial*> (matList->At(im));
-                 Int_t idmat = material->GetIndex();
-
-                 if (!global && idmat != proc->Medium()) continue;
-                 
-                 fMaterials[idmat] = im;
-                 // Skip media with no Cerenkov properties
-                 TFlukaCerenkov* cerenkovProp;
-                 if (!(cerenkovProp = dynamic_cast<TFlukaCerenkov*>(material->GetCerenkovProperties()))) continue;
-                 //
-                 // This medium has Cerenkov properties 
-                 //
-                 //
-                 // Write OPT-PROD card for each medium 
-                 Float_t  emin  = cerenkovProp->GetMinimumEnergy();
-                 Float_t  emax  = cerenkovProp->GetMaximumEnergy();          
-                 fprintf(pFlukaVmcInp, "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(pFlukaVmcInp, "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);
-                 
-                 if (cerenkovProp->IsMetal()) {
-                     fprintf(pFlukaVmcInp, "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(pFlukaVmcInp, "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(pFlukaVmcInp, "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(pFlukaVmcInp, "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);
-                 
-             } // materials
-         } else if (proc->Flag() == 0) {
-             fprintf(pFlukaVmcInp,"*\n*No Cerenkov photon generation\n");
-             fprintf(pFlukaVmcInp,"*Generated from call: SetProcess('CKOV',0)\n");
-             // zero = not used
-             // zero = not used
-             // zero = not used
-             // matMin = lower bound of the material indices in which the respective thresholds apply
-             // matMax = upper bound of the material indices in which the respective thresholds apply
-             // one = step length in assigning indices
-             //"CERE-OFF"; 
-             fprintf(pFlukaVmcInp,"OPT-PROD  %10.1f%10.1f%10.1f%10.1f%10.1f%10.1fCERE-OFF\n",zero,zero,zero,matMin,matMax,one);
-         }
-         else  {
-             fprintf(pFlukaVmcInp,"*\n*Illegal flag value in SetProcess('CKOV',?) call.\n");
-             fprintf(pFlukaVmcInp,"*No FLUKA card generated\n");
-         }
-      } // end of else if (strncmp(proc->GetName(),"CKOV",4) == 0)
-      
-      // Compton scattering
-      // G3 default value: 1
-      // G4 processes: G4ComptonScattering,
-      //               G4LowEnergyCompton,
-      //               G4PolarizedComptonScattering
-      // Particles: gamma
-      // Physics:   EM
-      // flag = 0 no Compton scattering
-      // flag = 1 Compton scattering, electron processed
-      // flag = 2 Compton scattering, no electron stored
-      // gMC ->SetProcess("COMP",1); // EMFCUT   -1.   0.  0. 3. lastmat 0. PHOT-THR
-      else if (strncmp(proc->GetName(),"COMP",4) == 0) {
-         if (proc->Flag() == 1 || proc->Flag() == 2) { 
-             fprintf(pFlukaVmcInp,"*\n*Energy threshold (GeV) for Compton scattering - resets to default=0.\n");
-             fprintf(pFlukaVmcInp,"*Generated from call: SetProcess('COMP',1);\n");
-             // - one = energy threshold (GeV) for Compton scattering - resets to default=0.
-             // zero = not used
-             // zero = not used
-             // matMin = lower bound of the material indices in which the respective thresholds apply
-             // matMax = upper bound of the material indices in which the respective thresholds apply
-             // one = step length in assigning indices
-             //"PHOT-THR"; 
-             theCut = -1.0;
-             nextc.Reset();
-             while ((cut = (TFlukaConfigOption*)nextc())) {
-                 if (strncmp(cut->GetName(), "CUTELE", 6) == 0 &&
-                     (cut->Medium() == proc->Medium())) theCut = cut->Cut();
-             }
-             fprintf(pFlukaVmcInp,"EMFCUT    %10.4g%10.1f%10.1f%10.1f%10.1f%10.1fPHOT-THR\n",theCut,zero,zero,matMin,matMax,one);
-         }
-         else if (proc->Flag() == 0) {
-             fprintf(pFlukaVmcInp,"*\n*No Compton scattering - no FLUKA card generated\n");
-             fprintf(pFlukaVmcInp,"*Generated from call: SetProcess('COMP',0)\n");
-         }
-         else  {
-             fprintf(pFlukaVmcInp,"*\n*Illegal flag value in SetProcess('COMP',?) call.\n");
-             fprintf(pFlukaVmcInp,"*No FLUKA card generated\n");
-         }
-      } // end of else if (strncmp(proc->GetName(),"COMP",4) == 0)
-      
-      // decay
-      // G3 default value: 1
-      // G4 process: G4Decay
-      // 
-      // Particles: all which decay is applicable for
-      // Physics:   General
-      // flag = 0 no decays
-      // flag = 1 decays, secondaries processed
-      // flag = 2 decays, no secondaries stored
-      //gMC ->SetProcess("DCAY",0); // not available
-      else if ((strncmp(proc->GetName(),"DCAY",4) == 0) && proc->Flag() == 0) 
-         cout << "SetProcess for flag =" << proc->GetName() << " value=" << proc->Flag() << " not avaliable!" << endl;
-      else if ((strncmp(proc->GetName(),"DCAY",4) == 0) && proc->Flag() == 1) {
-          // Nothing to do decays are switched on by default
-      }
-      
-      
-      // delta-ray
-      // G3 default value: 2
-      // !! G4 treats delta rays in different way
-      // G4 processes: G4eIonisation/G4IeIonization,
-      //               G4MuIonisation/G4IMuIonization,
-      //               G4hIonisation/G4IhIonisation
-      // Particles: charged
-      // Physics:   EM
-      // flag = 0 no energy loss
-      // flag = 1 restricted energy loss fluctuations
-      // flag = 2 complete energy loss fluctuations
-      // flag = 3 same as 1
-      // flag = 4 no energy loss fluctuations
-      // gMC ->SetProcess("DRAY",0); // DELTARAY 1.E+6 0.  0. 3. lastmat 0.
-      else if (strncmp(proc->GetName(),"DRAY",4) == 0) {
-         if (proc->Flag() == 0 || proc->Flag() == 4) {
-             fprintf(pFlukaVmcInp,"*\n*Kinetic energy threshold (GeV) for delta ray production\n");
-             fprintf(pFlukaVmcInp,"*Generated from call: SetProcess('DRAY',0) or SetProcess('DRAY',4)\n");
-             fprintf(pFlukaVmcInp,"*No delta ray production by muons - threshold set artificially high\n");
-             Double_t emin = 1.0e+6; // kinetic energy threshold (GeV) for delta ray production (discrete energy transfer)
-             // zero = ignored
-             // zero = ignored
-             // matMin = lower bound of the material indices in which the respective thresholds apply
-             // matMax = upper bound of the material indices in which the respective thresholds apply
-             // one = step length in assigning indices
-             fprintf(pFlukaVmcInp,"DELTARAY  %10.4g%10.1f%10.1f%10.1f%10.1f%10.1f\n",emin,zero,zero,matMin,matMax,one);
-         }
-         else if (proc->Flag() == 1 || proc->Flag() == 2 || proc->Flag() == 3) {
-             fprintf(pFlukaVmcInp,"*\n*Kinetic energy threshold (GeV) for delta ray production\n");
-             fprintf(pFlukaVmcInp,"*Generated from call: SetProcess('DRAY',flag), flag=1,2,3\n");
-             fprintf(pFlukaVmcInp,"*Delta ray production by muons switched on\n");
-             fprintf(pFlukaVmcInp,"*Energy threshold set by call SetCut('DCUTM',cut) or set to 1.0e+6.\n");
-             theCut = 1.0e+6;
-             nextc.Reset();
-             //
-             // Check cut one delta-rays from electrons
-             //
-             while ((cut = (TFlukaConfigOption*)nextc())) {
-                 if (strncmp(cut->GetName(), "DCUTM", 5) == 0 &&
-                     cut->Medium() == proc->Medium()) theCut = cut->Cut();
-             }
-             // theCut = kinetic energy threshold (GeV) for delta ray production (discrete energy transfer)
-             // zero = ignored
-             // zero = ignored
-             // matMin = lower bound of the material indices in which the respective thresholds apply
-             // matMax =  upper bound of the material indices in which the respective thresholds apply
-             // one = step length in assigning indices
-             fprintf(pFlukaVmcInp,"DELTARAY  %10.4g%10.1f%10.1f%10.1f%10.1f%10.1f\n",theCut,zero,zero,matMin,matMax,one);
-         }
-         else  {
-             fprintf(pFlukaVmcInp,"*\n*Illegal flag value in SetProcess('DRAY',?) call.\n");
-             fprintf(pFlukaVmcInp,"*No FLUKA card generated\n");
-         }
-      } // end of else if (strncmp(proc->GetName(),"DRAY",4) == 0)
-      
-      // hadronic process
-      // G3 default value: 1
-      // G4 processes: all defined by TG4PhysicsConstructorHadron
-      //  
-      // Particles: hadrons
-      // Physics:   Hadron
-      // flag = 0 no multiple scattering
-      // flag = 1 hadronic interactions, secondaries processed
-      // flag = 2 hadronic interactions, no secondaries stored
-      // gMC ->SetProcess("HADR",1); // ??? hadronic process
-      //Select pure GEANH (HADR 1) or GEANH/NUCRIN (HADR 3) ?????
-      else if (strncmp(proc->GetName(),"HADR",4) == 0) {
-         if (proc->Flag() == 1 || proc->Flag() == 2) {
-             fprintf(pFlukaVmcInp,"*\n*Hadronic interaction is ON by default in FLUKA\n");
-             fprintf(pFlukaVmcInp,"*No FLUKA card generated\n");
-         }
-         else if (proc->Flag() == 0) {
-             fprintf(pFlukaVmcInp,"*\n*Hadronic interaction is set OFF\n");
-             fprintf(pFlukaVmcInp,"*Generated from call: SetProcess('HADR',0);\n");
-             fprintf(pFlukaVmcInp,"*Switching off hadronic interactions not foreseen in FLUKA\n");
-             fprintf(pFlukaVmcInp,"THRESHOL  %10.1f%10.1f%10.1f%10.1e%10.1f\n",zero, zero, zero, 1.e10, zero);
-         }
-         else  {
-             fprintf(pFlukaVmcInp,"*\n*Illegal flag value in SetProcess('HADR',?) call.\n");
-             fprintf(pFlukaVmcInp,"*No FLUKA card generated\n");
-         }
-      } // end of else if (strncmp(proc->GetName(),"HADR",4) == 0)
-      
-      
-      // energy loss
-      // G3 default value: 2
-      // G4 processes: G4eIonisation/G4IeIonization,
-      //               G4MuIonisation/G4IMuIonization,
-      //               G4hIonisation/G4IhIonisation
-      // 
-      // Particles: charged
-      // Physics:   EM
-      // flag=0 no energy loss
-      // flag=1 restricted energy loss fluctuations
-      // flag=2 complete energy loss fluctuations
-      // flag=3 same as 1
-      // flag=4 no energy loss fluctuations
-      // If the value ILOSS is changed, then (in G3) cross-sections and energy
-      // loss tables must be recomputed via the command 'PHYSI'
-      // gMC ->SetProcess("LOSS",2); // ??? IONFLUCT ? energy loss
-      else if (strncmp(proc->GetName(),"LOSS",4) == 0) {
-         if (proc->Flag() > 0 || proc->Flag() < 4) { // restricted energy loss fluctuations
-             fprintf(pFlukaVmcInp,"*\n*Restricted energy loss fluctuations\n");
-             fprintf(pFlukaVmcInp,"*Generated from call: SetProcess('LOSS',1) or SetProcess('LOSS',3)\n");
-             // one = restricted energy loss fluctuations (for hadrons and muons) switched on
-             // one = restricted energy loss fluctuations (for e+ and e-) switched on
-             // one = minimal accuracy
-             // matMin = lower bound of the material indices in which the respective thresholds apply
-             // upper bound of the material indices in which the respective thresholds apply
-             fprintf(pFlukaVmcInp,"IONFLUCT  %10.1f%10.1f%10.1f%10.1f%10.1f\n",one, one, 4., matMin, matMax);
-         }
-         else if (proc->Flag() == 4) { // no energy loss fluctuations
-             fprintf(pFlukaVmcInp,"*\n*No energy loss fluctuations\n");
-             fprintf(pFlukaVmcInp,"*\n*Generated from call: SetProcess('LOSS',4)\n");
-             // - one = restricted energy loss fluctuations (for hadrons and muons) switched off
-             // - one = restricted energy loss fluctuations (for e+ and e-) switched off
-             // one = minimal accuracy
-             // matMin = lower bound of the material indices in which the respective thresholds apply
-             // matMax = upper bound of the material indices in which the respective thresholds apply
-             fprintf(pFlukaVmcInp,"IONFLUCT  %10.1f%10.1f%10.1f%10.1f%10.1f\n",-one,-one,one,matMin,matMax);
-         }
-         else  {
-             fprintf(pFlukaVmcInp,"*\n*Illegal flag value in SetProcess('LOSS',?) call.\n");
-             fprintf(pFlukaVmcInp,"*No FLUKA card generated\n");
-         }
-      } // end of else if (strncmp(proc->GetName(),"LOSS",4) == 0)
-      
-      
-      // multiple scattering
-      // G3 default value: 1
-      // G4 process: G4MultipleScattering/G4IMultipleScattering
-      // 
-      // Particles: charged
-      // Physics:   EM
-      // flag = 0 no multiple scattering
-      // flag = 1 Moliere or Coulomb scattering
-      // flag = 2 Moliere or Coulomb scattering
-      // flag = 3 Gaussian scattering
-      // gMC ->SetProcess("MULS",1); // MULSOPT multiple scattering
-      else if (strncmp(proc->GetName(),"MULS",4) == 0) {
-         if (proc->Flag() == 1 || proc->Flag() == 2 || proc->Flag() == 3) {
-             fprintf(pFlukaVmcInp,"*\n*Multiple scattering is ON by default for e+e- and for hadrons/muons\n");
-             fprintf(pFlukaVmcInp,"*No FLUKA card generated\n");
-         }
-         else if (proc->Flag() == 0) {
-             fprintf(pFlukaVmcInp,"*\n*Multiple scattering is set OFF\n");
-             fprintf(pFlukaVmcInp,"*Generated from call: SetProcess('MULS',0);\n");
-             // zero = ignored
-             // three = multiple scattering for hadrons and muons is completely suppressed
-             // three = multiple scattering for e+ and e- is completely suppressed
-             // matMin = lower bound of the material indices in which the respective thresholds apply
-             // matMax = upper bound of the material indices in which the respective thresholds apply
-             fprintf(pFlukaVmcInp,"MULSOPT   %10.1f%10.1f%10.1f%10.1f%10.1f\n",zero,three,three,matMin,matMax);
-         }
-         else  {
-             fprintf(pFlukaVmcInp,"*\n*Illegal flag value in SetProcess('MULS',?) call.\n");
-             fprintf(pFlukaVmcInp,"*No FLUKA card generated\n");
-         }
-      } // end of else if (strncmp(proc->GetName(),"MULS",4) == 0)
-      
-
-      // muon nuclear interaction
-      // G3 default value: 0
-      // G4 processes: G4MuNuclearInteraction,
-      // G4MuonMinusCaptureAtRest
-      // 
-      // Particles: mu
-      // Physics:   Not set
-      // flag = 0 no muon-nuclear interaction
-      // flag = 1 nuclear interaction, secondaries processed
-      // flag = 2 nuclear interaction, secondaries not processed
-      // gMC ->SetProcess("MUNU",1); // MUPHOTON  1.   0.  0. 3. lastmat
-      else if (strncmp(proc->GetName(),"MUNU",4) == 0) {
-         if (proc->Flag() == 1) {
-             fprintf(pFlukaVmcInp,"*\n*Muon nuclear interactions with production of secondary hadrons\n");
-             fprintf(pFlukaVmcInp,"*\n*Generated from call: SetProcess('MUNU',1);\n");
-             // one = full simulation of muon nuclear interactions and production of secondary hadrons
-             // zero = ratio of longitudinal to transverse virtual photon cross-section - Default = 0.25.
-             // zero = fraction of rho-like interactions ( must be < 1) - Default = 0.75.
-             // matMin = lower bound of the material indices in which the respective thresholds apply
-             // matMax = upper bound of the material indices in which the respective thresholds apply
-             fprintf(pFlukaVmcInp,"MUPHOTON  %10.1f%10.1f%10.1f%10.1f%10.1f\n",one,zero,zero,matMin,matMax);
-         }
-         else if (proc->Flag() == 2) {
-             fprintf(pFlukaVmcInp,"*\n*Muon nuclear interactions without production of secondary hadrons\n");
-             fprintf(pFlukaVmcInp,"*Generated from call: SetProcess('MUNU',2);\n");
-             // two = full simulation of muon nuclear interactions and production of secondary hadrons
-             // zero = ratio of longitudinal to transverse virtual photon cross-section - Default = 0.25.
-             // zero = fraction of rho-like interactions ( must be < 1) - Default = 0.75.
-             // matMin = lower bound of the material indices in which the respective thresholds apply
-             // matMax = upper bound of the material indices in which the respective thresholds apply
-             fprintf(pFlukaVmcInp,"MUPHOTON  %10.1f%10.1f%10.1f%10.1f%10.1f\n",two,zero,zero,matMin,matMax);
-         }
-         else if (proc->Flag() == 0) {
-             fprintf(pFlukaVmcInp,"*\n*No muon nuclear interaction - no FLUKA card generated\n");
-             fprintf(pFlukaVmcInp,"*Generated from call: SetProcess('MUNU',0)\n");
-         }
-         else  {
-             fprintf(pFlukaVmcInp,"*\n*Illegal flag value in SetProcess('MUNU',?) call.\n");
-             fprintf(pFlukaVmcInp,"*No FLUKA card generated\n");
-         }
-      } // end of else if (strncmp(proc->GetName(),"MUNU",4) == 0)
-      
-      
-      // photofission
-      // G3 default value: 0
-      // G4 process: ??
-      //
-      // Particles: gamma
-      // Physics:   ??
-      // gMC ->SetProcess("PFIS",0); // PHOTONUC -1.   0.  0. 3. lastmat 0.
-      // flag = 0 no photon fission
-      // flag = 1 photon fission, secondaries processed
-      // flag = 2 photon fission, no secondaries stored
-      else if (strncmp(proc->GetName(),"PFIS",4) == 0) {
-         if (proc->Flag() == 0) {
-             fprintf(pFlukaVmcInp,"*\n*No photonuclear interactions\n");
-             fprintf(pFlukaVmcInp,"*Generated from call: SetProcess('PFIS',0);\n");
-             // - one = no photonuclear interactions
-             // zero = not used
-             // zero = not used
-             // matMin = lower bound of the material indices in which the respective thresholds apply
-             // matMax = upper bound of the material indices in which the respective thresholds apply
-             fprintf(pFlukaVmcInp,"PHOTONUC  %10.1f%10.1f%10.1f%10.1f%10.1f\n",-one,zero,zero,matMin,matMax);
-         }
-         else if (proc->Flag() == 1) {
-             fprintf(pFlukaVmcInp,"*\n*Photon nuclear interactions are activated at all energies\n");
-             fprintf(pFlukaVmcInp,"*Generated from call: SetProcess('PFIS',1);\n");
-             // one = photonuclear interactions are activated at all energies
-             // zero = not used
-             // zero = not used
-             // matMin = lower bound of the material indices in which the respective thresholds apply
-             // matMax = upper bound of the material indices in which the respective thresholds apply
-             fprintf(pFlukaVmcInp,"PHOTONUC  %10.1f%10.1f%10.1f%10.1f%10.1f\n",one,zero,zero,matMin,matMax);
-         }
-         else if (proc->Flag() == 0) {
-             fprintf(pFlukaVmcInp,"*\n*No photofission - no FLUKA card generated\n");
-             fprintf(pFlukaVmcInp,"*Generated from call: SetProcess('PFIS',0)\n");
-         }
-         else {
-             fprintf(pFlukaVmcInp,"*\n*Illegal flag value in SetProcess('PFIS',?) call.\n");
-             fprintf(pFlukaVmcInp,"*No FLUKA card generated\n");
-         }
-      }
-
-      // photo electric effect
-      // G3 default value: 1
-      // G4 processes: G4PhotoElectricEffect
-      //               G4LowEnergyPhotoElectric
-      // Particles: gamma
-      // Physics:   EM
-      // flag = 0 no photo electric effect
-      // flag = 1 photo electric effect, electron processed
-      // flag = 2 photo electric effect, no electron stored
-      // gMC ->SetProcess("PHOT",1); // EMFCUT    0.  -1.  0. 3. lastmat 0. PHOT-THR
-      else if (strncmp(proc->GetName(),"PHOT",4) == 0) {
-         if (proc->Flag() == 1) {
-             fprintf(pFlukaVmcInp,"*\n*Photo electric effect is activated\n");
-             fprintf(pFlukaVmcInp,"*Generated from call: SetProcess('PHOT',1);\n");
-             // zero = ignored
-             // - one = resets to default=0.
-             // zero = ignored
-             // matMin = lower bound of the material indices in which the respective thresholds apply
-             // matMax = upper bound of the material indices in which the respective thresholds apply
-             // one = step length in assigning indices
-             //"PHOT-THR"; 
-             theCut = -1.0;
-             nextc.Reset();
-             while ((cut = (TFlukaConfigOption*)nextc())) {
-                 if (strncmp(cut->GetName(), "CUTELE", 6) == 0 &&
-                     (cut->Medium() == proc->Medium())) theCut = cut->Cut();
-             }
-             fprintf(pFlukaVmcInp,"EMFCUT    %10.1f%10.4g%10.1f%10.1f%10.1f%10.1fPHOT-THR\n",zero,theCut,zero,matMin,matMax,one);
-         }
-         else if (proc->Flag() == 0) {
-             fprintf(pFlukaVmcInp,"*\n*No photo electric effect - no FLUKA card generated\n");
-             fprintf(pFlukaVmcInp,"*Generated from call: SetProcess('PHOT',0)\n");
-         }
-         else {
-             fprintf(pFlukaVmcInp,"*\n*Illegal flag value in SetProcess('PHOT',?) call.\n");
-             fprintf(pFlukaVmcInp,"*No FLUKA card generated\n");
-         }
-      } // else if (strncmp(proc->GetName(),"PHOT",4) == 0)
-      
-      
-      // Rayleigh scattering
-      // G3 default value: 0
-      // G4 process: G4OpRayleigh
-      // 
-      // Particles: optical photon
-      // Physics:   Optical
-      // flag = 0 Rayleigh scattering off
-      // flag = 1 Rayleigh scattering on
-      //xx gMC ->SetProcess("RAYL",1);
-      else if (strncmp(proc->GetName(),"RAYL",4) == 0) {
-         if (proc->Flag() == 1) {
-             fprintf(pFlukaVmcInp,"*\n*Rayleigh scattering is ON by default in FLUKA\n");
-             fprintf(pFlukaVmcInp,"*No FLUKA card generated\n");
-         }
-         else if (proc->Flag() == 0) {
-             fprintf(pFlukaVmcInp,"*\n*Rayleigh scattering is set OFF\n");
-             fprintf(pFlukaVmcInp,"*Generated from call: SetProcess('RAYL',0);\n");
-             // - one = no Rayleigh scattering and no binding corrections for Compton
-             // matMin = lower bound of the material indices in which the respective thresholds apply
-             // matMax = upper bound of the material indices in which the respective thresholds apply
-             fprintf(pFlukaVmcInp,"EMFRAY    %10.1f%10.1f%10.1f%10.1f\n",-one,three,matMin,matMax);
-         }
-         else  {
-             fprintf(pFlukaVmcInp,"*\n*Illegal flag value in SetProcess('RAYL',?) call.\n");
-             fprintf(pFlukaVmcInp,"*No FLUKA card generated\n");
-         }
-      } // end of else if (strncmp(proc->GetName(),"RAYL",4) == 0)
-      
-      
-      // 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(proc->GetName(),"SYNC",4) == 0) {
-         fprintf(pFlukaVmcInp,"*\n*Synchrotron radiation generation is NOT implemented in FLUKA\n");
-         fprintf(pFlukaVmcInp,"*No FLUKA card generated\n");
-      }
-      
-      
-      // 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(proc->GetName(),"AUTO",4) == 0) {
-         fprintf(pFlukaVmcInp,"*\n*Automatic calculation of tracking medium parameters is always ON in FLUKA\n");
-         fprintf(pFlukaVmcInp,"*No FLUKA card generated\n");
-      }
-      
-      
-      // 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(proc->GetName(),"STRA",4) == 0) {
-         if (proc->Flag() == 0 || proc->Flag() == 2 || proc->Flag() == 3) {
-             fprintf(pFlukaVmcInp,"*\n*Ionization energy losses calculation is activated\n");
-             fprintf(pFlukaVmcInp,"*Generated from call: SetProcess('STRA',n);, n=0,1,2\n");
-             // one = restricted energy loss fluctuations (for hadrons and muons) switched on
-             // one = restricted energy loss fluctuations (for e+ and e-) switched on
-             // one = minimal accuracy
-             // matMin = lower bound of the material indices in which the respective thresholds apply
-             // matMax = upper bound of the material indices in which the respective thresholds apply
-             fprintf(pFlukaVmcInp,"IONFLUCT  %10.1f%10.1f%10.1f%10.1f%10.1f\n",one,one,one,matMin,matMax);
-         }
-         else {
-             fprintf(pFlukaVmcInp,"*\n*Illegal flag value in SetProcess('STRA',?) call.\n");
-             fprintf(pFlukaVmcInp,"*No FLUKA card generated\n");
-         }
-      } // else if (strncmp(proc->GetName(),"STRA",4) == 0)
-      
-
-
-
-      else { // processes not yet treated
-         
-         // light photon absorption (Cerenkov photons)
-         // it is turned on when Cerenkov process is turned on
-         // G3 default value: 0
-         // G4 process: G4OpAbsorption, G4OpBoundaryProcess
-         // 
-         // Particles: optical photon
-         // Physics:   Optical
-         // flag = 0 no absorption of Cerenkov photons
-         // flag = 1 absorption of Cerenkov photons
-         // gMC ->SetProcess("LABS",2); // ??? Cerenkov light absorption
-         
-
-
-         cout << "SetProcess for flag=" << proc->GetName() << " value=" << proc->Flag() << " not yet implemented!" << endl;
-      }
-  } //end of loop number of SetProcess calls
-
-// Loop over number of SetCut calls  
-           
-  nextc.Reset();
-  while ((cut = (TFlukaConfigOption*)nextc())) {
-      Float_t matMin = three;
-      Float_t matMax = fLastMaterial;
-      Bool_t global  = kTRUE;
-      if (cut->Medium() != -1) {
-         Int_t mat;
-         if ((mat = cut->Medium()) >= GetFlukaMaterials()->GetEntries()) continue;
-         matMin = Float_t(mat);
-         matMax = matMin;
-         global = kFALSE;
-         TGeoMaterial*    material =  (TGeoMaterial*) (GetFlukaMaterials())->At(GetMaterialIndex(mat));
-         fprintf(pFlukaVmcInp,"*\n*Material specific cut setting for #%8d %s %s %13.3e\n", 
-                 mat, material->GetName(), cut->GetName(), cut->Cut());  
-         
-      } 
-      
-      // cuts handled in SetProcess calls
-      if      (strncmp(cut->GetName(),"BCUTM",5)  == 0) continue;
-      else if (strncmp(cut->GetName(),"BCUTE",5)  == 0) continue;
-      else if (strncmp(cut->GetName(),"DCUTM",5)  == 0) continue;
-      else if (strncmp(cut->GetName(),"PPCUTM",6) == 0) continue;
-      
-      // gammas
-      // G4 particles: "gamma"
-      // G3 default value: 0.001 GeV
-      // gMC ->SetCut("CUTGAM",cut); // cut for gammas
-      
-      else if (strncmp(cut->GetName(),"CUTGAM",6) == 0 && global) {
-         fprintf(pFlukaVmcInp,"*\n*Cut for gamma\n");
-         fprintf(pFlukaVmcInp,"*Generated from call: SetCut('CUTGAM',cut);\n");
-         fprintf(pFlukaVmcInp,"EMFCUT    %10.4g%10.4g%10.1f%10.1f%10.1f%10.1f\n", 
-                 zero, cut->Cut(), zero, zero, Float_t(fGeom->NofVolumes()), one);
-      }
-      else if (strncmp(cut->GetName(),"CUTGAM",6) == 0 && !global) {
-         // loop over materials for EMFCUT FLUKA cards
-         for (j=0; j < matMax-matMin+1; j++) {
-             Int_t nreg, imat, *reglist;
-             Float_t ireg;
-             imat = (Int_t) matMin + j;
-             reglist = fGeom->GetMaterialList(imat, nreg);
-             // loop over regions of a given material
-             for (Int_t k = 0; k < nreg; k++) {
-                 ireg = reglist[k];
-                 fprintf(pFlukaVmcInp,"EMFCUT    %10.4g%10.4g%10.1f%10.1f%10.1f%10.1f\n", zero, cut->Cut(), zero, ireg, ireg, one);
-             }
-         }
-      } // end of else if for gamma
-      
-
-      // electrons
-      // G4 particles: "e-"
-      // ?? positrons
-      // G3 default value: 0.001 GeV
-      //gMC ->SetCut("CUTELE",cut); // cut for e+,e-
-      else if (strncmp(cut->GetName(),"CUTELE",6) == 0 && global) {
-         fprintf(pFlukaVmcInp,"*\n*Cut for electrons\n");
-         fprintf(pFlukaVmcInp,"*Generated from call: SetCut('CUTELE',cut);\n");
-         fprintf(pFlukaVmcInp,"EMFCUT    %10.4g%10.4g%10.1f%10.1f%10.1f%10.1f\n", 
-                 -cut->Cut(), zero, zero, zero, Float_t(fGeom->NofVolumes()), one);
-      }
-      else if (strncmp(cut->GetName(),"CUTELE",6) == 0 && !global) {
-         // loop over materials for EMFCUT FLUKA cards
-         for (j=0; j < matMax-matMin+1; j++) {
-             Int_t nreg, imat, *reglist;
-             Float_t ireg;
-             imat = (Int_t) matMin + j;
-             reglist = fGeom->GetMaterialList(imat, nreg);
-             // loop over regions of a given material
-             for (k=0; k<nreg; k++) {
-                 ireg = reglist[k];
-                 fprintf(pFlukaVmcInp,"EMFCUT    %10.4g%10.4g%10.1f%10.1f%10.1f%10.1f\n", -cut->Cut(), zero, zero, ireg, ireg, one);
-             }
-         }
-      } // end of else if for electrons
+    } //end of while until START card
+    
+ fin:
 
+    // Pass information to configuration objects
     
-      // neutral hadrons
-      // G4 particles: of type "baryon", "meson", "nucleus" with zero charge
-      // G3 default value: 0.01 GeV
-      //gMC ->SetCut("CUTNEU",cut); // cut for neutral hadrons
-      else if (strncmp(cut->GetName(),"CUTNEU",6) == 0 && global) {
-       fprintf(pFlukaVmcInp,"*\n*Cut for neutral hadrons\n");
-       fprintf(pFlukaVmcInp,"*Generated from call: SetCut('CUTNEU',cut);\n");
-         
-       // 8.0 = Neutron
-       // 9.0 = Antineutron
-       fprintf(pFlukaVmcInp,"PART-THR  %10.4g%10.1f%10.1f\n",-cut->Cut(),8.0,9.0);
-         
-       // 12.0 = Kaon zero long
-       // 12.0 = Kaon zero long
-       fprintf(pFlukaVmcInp,"PART-THR  %10.4g%10.1f%10.1f\n",-cut->Cut(),12.0,12.0);
-         
-       // 17.0 = Lambda, 18.0 = Antilambda
-       // 19.0 = Kaon zero short
-       fprintf(pFlukaVmcInp,"PART-THR  %10.4g%10.1f%10.1f\n",-cut->Cut(),17.0,19.0);
-         
-       // 22.0 = Sigma zero, Pion zero, Kaon zero
-       // 25.0 = Antikaon zero
-       fprintf(pFlukaVmcInp,"PART-THR  %10.4g%10.1f%10.1f\n",-cut->Cut(),22.0,25.0);
-         
-       // 32.0 = Antisigma zero
-       // 32.0 = Antisigma zero
-       fprintf(pFlukaVmcInp,"PART-THR  %10.4g%10.1f%10.1f\n",-cut->Cut(),32.0,32.0);
-         
-       // 34.0 = Xi zero
-       // 35.0 = AntiXi zero
-       fprintf(pFlukaVmcInp,"PART-THR  %10.4g%10.1f%10.1f\n",-cut->Cut(),34.0,35.0);
-         
-       // 47.0 = D zero
-       // 48.0 = AntiD zero
-       fprintf(pFlukaVmcInp,"PART-THR  %10.4g%10.1f%10.1f\n",-cut->Cut(),47.0,48.0);
-         
-       // 53.0 = Xi_c zero
-       // 53.0 = Xi_c zero
-       fprintf(pFlukaVmcInp,"PART-THR  %10.4g%10.1f%10.1f\n",-cut->Cut(),53.0,53.0);
-         
-       // 55.0 = Xi'_c zero
-       // 56.0 = Omega_c zero
-       fprintf(pFlukaVmcInp,"PART-THR  %10.4g%10.1f%10.1f\n",-cut->Cut(),55.0,56.0);
-         
-       // 59.0 = AntiXi_c zero
-       // 59.0 = AntiXi_c zero
-       fprintf(pFlukaVmcInp,"PART-THR  %10.4g%10.1f%10.1f\n",-cut->Cut(),59.0,59.0);
-         
-       // 61.0 = AntiXi'_c zero
-       // 62.0 = AntiOmega_c zero
-       fprintf(pFlukaVmcInp,"PART-THR  %10.4g%10.1f%10.1f\n",-cut->Cut(),61.0,62.0);
-      }
-      
-      // charged hadrons
-      // G4 particles: of type "baryon", "meson", "nucleus" with non-zero charge
-      // G3 default value: 0.01 GeV
-      //gMC ->SetCut("CUTHAD",cut); // cut for charged hadrons
-      else if (strncmp(cut->GetName(),"CUTHAD",6) == 0 && global) {
-       fprintf(pFlukaVmcInp,"*\n*Cut for charged hadrons\n");
-       fprintf(pFlukaVmcInp,"*Generated from call: SetCut('CUTHAD',cut);\n");
-         
-       // 1.0 = Proton
-       // 2.0 = Antiproton
-       fprintf(pFlukaVmcInp,"PART-THR  %10.4g%10.1f%10.1f\n",-cut->Cut(),1.0,2.0);
-         
-       // 13.0 = Positive Pion, Negative Pion, Positive Kaon
-       // 16.0 = Negative Kaon
-       fprintf(pFlukaVmcInp,"PART-THR  %10.4g%10.1f%10.1f\n",-cut->Cut(),13.0,16.0);
-         
-       // 20.0 = Negative Sigma
-       // 21.0 = Positive Sigma
-       fprintf(pFlukaVmcInp,"PART-THR  %10.4g%10.1f%10.1f\n",-cut->Cut(),20.0,21.0);
-         
-       // 31.0 = Antisigma minus
-       // 33.0 = Antisigma plus
-       // 2.0 = step length
-       fprintf(pFlukaVmcInp,"PART-THR  %10.4g%10.1f%10.1f%10.1f\n",-cut->Cut(),31.0,33.0,2.0);
-         
-       // 36.0 = Negative Xi, Positive Xi, Omega minus
-       // 39.0 = Antiomega
-       fprintf(pFlukaVmcInp,"PART-THR  %10.4g%10.1f%10.1f\n",-cut->Cut(),36.0,39.0);
-         
-       // 45.0 = D plus
-       // 46.0 = D minus
-       fprintf(pFlukaVmcInp,"PART-THR  %10.4g%10.1f%10.1f\n",-cut->Cut(),45.0,46.0);
-         
-       // 49.0 = D_s plus, D_s minus, Lambda_c plus
-       // 52.0 = Xi_c plus
-       fprintf(pFlukaVmcInp,"PART-THR  %10.4g%10.1f%10.1f\n",-cut->Cut(),49.0,52.0);
-         
-       // 54.0 = Xi'_c plus
-       // 60.0 = AntiXi'_c minus
-       // 6.0 = step length
-       fprintf(pFlukaVmcInp,"PART-THR  %10.4g%10.1f%10.1f%10.1f\n",-cut->Cut(),54.0,60.0,6.0);
-         
-       // 57.0 = Antilambda_c minus
-       // 58.0 = AntiXi_c minus
-       fprintf(pFlukaVmcInp,"PART-THR  %10.4g%10.1f%10.1f\n",-cut->Cut(),57.0,58.0);
-      }
-
-      // muons
-      // G4 particles: "mu+", "mu-"
-      // G3 default value: 0.01 GeV
-      //gMC ->SetCut("CUTMUO",cut); // cut for mu+, mu-
-      else if (strncmp(cut->GetName(),"CUTMUO",6)== 0 && global) {
-       fprintf(pFlukaVmcInp,"*\n*Cut for muons\n");
-       fprintf(pFlukaVmcInp,"*Generated from call: SetCut('CUTMUO',cut);\n");
-       // 10.0 = Muon+
-       // 11.0 = Muon-
-       fprintf(pFlukaVmcInp,"PART-THR  %10.4g%10.1f%10.1f\n",-cut->Cut(),10.0,11.0);
-      }
-      
-      //
-      // time of flight cut in seconds
-      // G4 particles: all
-      // G3 default value: 0.01 GeV
-      //gMC ->SetCut("TOFMAX",tofmax); // time of flight cuts in seconds
-      else if (strncmp(cut->GetName(),"TOFMAX",6) == 0) {
-       fprintf(pFlukaVmcInp,"*\n*Time of flight cuts in seconds\n");
-       fprintf(pFlukaVmcInp,"*Generated from call: SetCut('TOFMAX',tofmax);\n");
-       // zero = ignored
-       // zero = ignored
-       // -6.0 = lower bound of the particle numbers for which the transport time cut-off and/or the start signal is to be applied
-       // 64.0 = upper bound of the particle numbers for which the transport time cut-off and/or the start signal is to be applied
-       fprintf(pFlukaVmcInp,"TIME-CUT  %10.4g%10.1f%10.1f%10.1f%10.1f\n",cut->Cut()*1.e9,zero,zero,-6.0,64.0);
-      }
-      
-      else if (global){
-       cout << "SetCut for flag=" << cut->GetName() << " value=" << cut->Cut() << " not yet implemented!" << endl;
-      }
-      else {
-       cout << "SetCut for flag=" << cut->GetName() << " value=" << cut->Cut() << " (material specific) not yet implemented!" << endl;
-      }
-      
-  } //end of loop over SetCut calls
-  
+    Float_t fLastMaterial = fGeom->GetLastMaterialIndex();
+    TFlukaConfigOption::SetStaticInfo(pFlukaVmcInp, 3, fLastMaterial, fGeom);
+    
+    TIter next(fUserConfig);
+    TFlukaConfigOption* proc;
+    while((proc = (TFlukaConfigOption*)next())) proc->WriteFlukaInputCards();
+
 // Add START and STOP card
   fprintf(pFlukaVmcInp,"START     %10.1f\n",fEventsPerRun);
   fprintf(pFlukaVmcInp,"STOP      \n");
    
   
 // Close files
-  
    fclose(pFlukaVmcCoreInp);
    fclose(pFlukaVmcFlukaMat);
    fclose(pFlukaVmcInp);
-   
+
+
+//
+// Initialisation needed for Cerenkov photon production and transport
+    TObjArray *matList = GetFlukaMaterials();
+    Int_t nmaterial =  matList->GetEntriesFast();
+    fMaterials = new Int_t[nmaterial+3];
+    
+    for (Int_t im = 0; im < nmaterial; im++)
+    {
+       TGeoMaterial* material = dynamic_cast<TGeoMaterial*> (matList->At(im));
+       Int_t idmat = material->GetIndex();
+       fMaterials[idmat] = im;
+    }
 } // end of InitPhysics
 
 
index 2a117102635b41f4105da84e1e46aa034c9863d3..f40cee50542d9822e8f9f5baff31e63596410aab 100644 (file)
@@ -123,8 +123,7 @@ class TFluka : public TVirtualMC {
   virtual void       SetProcess(const char* flagName, Int_t flagValue, Int_t imed);
   virtual Bool_t     SetCut(const char* cutName, Double_t cutValue);
   virtual void       SetCut(const char* cutName, Double_t cutValue, Int_t imed);
-  virtual TObjArray* GetListOfProceses() {return fProcesses;}
-  virtual TObjArray* GetListOfCuts()     {return fCuts;}
+  virtual TObjArray* GetListOfUserConfigs() {return fUserConfig;}
   virtual Double_t   Xsec(char*, Double_t, Int_t, Int_t);
   
   // Particle table usage         
@@ -377,8 +376,7 @@ class TFluka : public TVirtualMC {
   TGeoMCGeometry      *fMCGeo;              // Interface to TGeo builder
 
   // SetProcess, SetCut and user Scoring dynamic storage
-  TObjArray* fProcesses;             // List of processes
-  TObjArray* fCuts;                  // List of cuts
+  TObjArray* fUserConfig;            // List of user physics configuration 
   TObjArray* fUserScore;             // List of user scoring options
   
   ClassDef(TFluka,1)  //C++ interface to Fluka montecarlo
index d2eecff94f6aeb1df362602c06dff73acbecbce8..7c39dff54f668fd31818b112afd7b7cfe7ef3a7f 100644 (file)
 /* $Id$*/
 
 #include "TFlukaConfigOption.h"
+#include "TFlukaMCGeometry.h"
+#include "TFluka.h"
+#include "TFlukaCerenkov.h"
+
+#include <TString.h>
+#include <TObjArray.h>
+#include <TVirtualMC.h>
+#include <TGeoMaterial.h>
+
+Float_t            TFlukaConfigOption::fgMatMin(-1.);
+Float_t            TFlukaConfigOption::fgMatMax(-1.);
+FILE*              TFlukaConfigOption::fgFile(0x0);
+TFlukaMCGeometry*  TFlukaConfigOption::fgGeom(0x0);
+
+Double_t TFlukaConfigOption::fgDCutValue[11];
+Int_t    TFlukaConfigOption::fgDProcessFlag[15];
+
+
 ClassImp(TFlukaConfigOption)
 
 
 TFlukaConfigOption::TFlukaConfigOption()
 {
     // Default constructor
+    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;
 }
 
 
-TFlukaConfigOption::TFlukaConfigOption(const char* cutName, Double_t cut)
-    : TNamed(cutName, "Cut")
+TFlukaConfigOption::TFlukaConfigOption(Int_t medium)
 {
     // Constructor
-    fCutValue        = cut;
-    fMedium          = -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;
 }
 
-TFlukaConfigOption::TFlukaConfigOption(const char* cutName, Double_t cut, Int_t imed)
-    : TNamed(cutName, "Cut")
+void TFlukaConfigOption::SetCut(const char* flagname, Double_t val)
 {
-    // Constructor
-    fCutValue       = cut;
-    fMedium         = imed;
+    // Set a cut value
+    const TString cuts[11] = 
+       {"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;
+       }
+    }
 }
 
+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"};
+    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;
+       }
+    }
+}
 
-TFlukaConfigOption::TFlukaConfigOption(const char* procName, Int_t flag)
-    : TNamed(procName, "Process")
+void TFlukaConfigOption::WriteFlukaInputCards()
 {
-    // Constructor
-    fProcessFlag     = flag;
-    fMedium          = -1;
+    // 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;
+    } else {
+       fprintf(fgFile,"*\n*Global process and cut settings \n");
+       fCMatMin = fgMatMin;
+       fCMatMax = fgMatMax;
+    }
+
+//
+// Handle Process Flags 
+//    
+    if (fProcessFlag[kDCAY] != -1) ProcessDCAY();
+    if (fProcessFlag[kPAIR] != -1) ProcessPAIR();
+    if (fProcessFlag[kBREM] != -1) ProcessBREM();
+    if (fProcessFlag[kCOMP] != -1) ProcessCOMP();
+    if (fProcessFlag[kPHOT] != -1) ProcessPHOT();
+    if (fProcessFlag[kPFIS] != -1) ProcessPFIS();
+    if (fProcessFlag[kANNI] != -1) ProcessANNI();
+    if (fProcessFlag[kMUNU] != -1) ProcessMUNU();
+    if (fProcessFlag[kHADR] != -1) ProcessHADR();
+    if (fProcessFlag[kMULS] != -1) ProcessMULS();
+    if (fProcessFlag[kRAYL] != -1) ProcessRAYL();
+
+    if (fProcessFlag[kLOSS] != -1 || fProcessFlag[kDRAY] != -1)                                    ProcessLOSS();
+    if ((fMedium == -1 && fProcessFlag[kCKOV] > 0) || (fMedium > -1 && fProcessFlag[kCKOV] != -1)) ProcessCKOV();
+    
+//
+// Handle Cuts
+//
+    if (fCutValue[kCUTGAM] >= 0.) ProcessCUTGAM();
+    if (fCutValue[kCUTELE] >= 0.) ProcessCUTELE();
+    if (fCutValue[kCUTNEU] >= 0.) ProcessCUTNEU();
+    if (fCutValue[kCUTHAD] >= 0.) ProcessCUTHAD();
+    if (fCutValue[kCUTMUO] >= 0.) ProcessCUTMUO();
+
+    if (fCutValue[kTOFMAX] >= 0.) ProcessTOFMAX();
 }
 
-TFlukaConfigOption::TFlukaConfigOption(const char* procName, Int_t flag, Int_t imed)
-    : TNamed(procName, "Process")
+void TFlukaConfigOption::ProcessDCAY()
 {
-    // Constructor
-    fProcessFlag    = flag;
-    fMedium         = imed;
+    // Process DCAY option
+    fprintf(fgFile,"*\n* --- DCAY --- Decays. Flag = %5d\n", fProcessFlag[kDCAY]);
+    if (fProcessFlag[kDCAY] == 0) {
+       printf("Decays cannot be switched off \n");
+    } else {
+       fprintf(fgFile, "* Decays are on by default\n");
+    }
+}
+
+
+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]);
+    //
+    // 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.);
+    } else {
+       fprintf(fgFile,"EMFCUT    %10.1f%10.1f%10.4g%10.1f%10.1f%10.1fPHOT-THR\n",0., 0., 1e10,  fCMatMin, fCMatMax, 1.);
+    }
+    
+    //
+    // Direct pair production by Muons and Hadrons
+    //
+    //
+    // Attention ! This card interferes with BREM
+    //
+
+    if (fProcessFlag[kBREM] == -1 ) fProcessFlag[kBREM] = fgDProcessFlag[kBREM];
+    if (fCutValue[kBCUTM]   == -1.) fCutValue[kBCUTM]   = fgDCutValue[kBCUTM]; 
+
+
+    Float_t flag = -3.;    
+    if (fProcessFlag[kPAIR] >  0 && fProcessFlag[kBREM] == 0) flag =  1.;
+    if (fProcessFlag[kPAIR] == 0 && fProcessFlag[kBREM]  > 0) flag =  2.;
+    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;
+    
+    //
+    // Energy cut for pair prodution
+    //
+    Float_t cutP = fCutValue[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;
+    if (cutP < 0.) cutP = 0.;
+    // No explicite generation of e+/e-
+    if (fProcessFlag[kPAIR] == 2) cutP = -1.;
+    //
+    // Energy cut for bremsstrahlung
+    //
+    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]);
+
+       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);
+}
+
+
+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]);
+
+    //
+    // e+/- -> e+/- gamma
+    //
+    Float_t cutB = fCutValue[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.);
+    } 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.);
+    }
+
+    //
+    // Bremsstrahlung by muons and hadrons
+    //
+    cutB = fCutValue[kBCUTM];
+    if (fCutValue[kBCUTM]   == -1.) cutB = fgDCutValue[kBCUTM];        
+    if (fProcessFlag[kBREM] == 2) cutB = -1.;
+    Float_t flag = 2.;
+    if (fProcessFlag[kBREM] == 0) flag = -2.;
+    
+    fprintf(fgFile,"PAIRBREM  %10.1f%10.4g%10.4g%10.1f%10.1f\n", flag, 0., cutB, fCMatMin, fCMatMax);
+}
+
+void TFlukaConfigOption::ProcessCOMP()
+{
+    // Process COMP option
+    fprintf(fgFile,"*\n* --- COMP --- Compton scattering  Flag = %5d \n", fProcessFlag[kCOMP]);
+
+    //
+    // Compton scattering
+    //
+
+    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.);
+    } 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.);
+    }
+}
+
+void TFlukaConfigOption::ProcessPHOT()
+{
+    // Process PHOS option
+    fprintf(fgFile,"*\n* --- PHOT --- Photoelectric effect. Flag = %5d\n", fProcessFlag[kPHOT]);
+
+    //
+    // Photoelectric effect
+    //
+
+    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.);
+    } 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.);
+    }
+}
+
+void TFlukaConfigOption::ProcessANNI()
+{
+    // Process ANNI option
+    fprintf(fgFile,"*\n* --- ANNI --- Positron annihilation. Flag = %5d \n", fProcessFlag[kANNI]);
+
+    //
+    // Positron annihilation
+    //
+
+    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.);
+    } 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.);
+    }
+}
+
+
+void TFlukaConfigOption::ProcessPFIS()
+{
+    // Process PFIS option
+    fprintf(fgFile,"*\n* --- PFIS --- Photonuclear interaction  Flag = %5d\n", fProcessFlag[kPFIS]);
+
+    //
+    // Photonuclear interactions
+    //
+
+    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.);
+    } else {
+       fprintf(fgFile,"PHOTONUC  %10.1f%10.1f%10.1f%10.1f%10.1f%10.1f\n",-1.                          , 0., 0.,  fCMatMin, fCMatMax, 1.);
+    }
+}
+
+void TFlukaConfigOption::ProcessMUNU()
+{
+    // Process MUNU option
+    fprintf(fgFile,"*\n* --- MUNU --- Muon nuclear interaction. Flag = %5d\n", fProcessFlag[kMUNU]);
+    
+    //
+    // 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.);
+    } else {
+       fprintf(fgFile,"MUPHOTON  %10.1f%10.1f%10.1f%10.1f%10.1f%10.1f\n",-1.                          , 0.,   0.,    fCMatMin, fCMatMax, 1.);
+    }
+}
+
+void TFlukaConfigOption::ProcessRAYL()
+{
+    // Process RAYL option
+    fprintf(fgFile,"*\n* --- RAYL --- Rayleigh Scattering. Flag = %5d\n", fProcessFlag[kRAYL]);
+
+    //
+    // Rayleigh scattering
+    //
+    Int_t nreg;
+    Int_t* reglist = fgGeom->GetMaterialList(fMedium, nreg);
+    //
+    // 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.);
+       }
+    }
+}
+
+void TFlukaConfigOption::ProcessCKOV()
+{
+    // Process CKOV option
+    fprintf(fgFile,"*\n* --- CKOV --- Cerenkov Photon production.  %5d\n", fProcessFlag[kCKOV]);
+
+    //
+    // Cerenkov photon production
+    //
+
+    TFluka* fluka = (TFluka*) gMC;
+    TObjArray *matList = fluka->GetFlukaMaterials();
+    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();
+//
+// 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);
+           
+           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-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);
+       } 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.);
+       }
+    }
+}
+
+
+void TFlukaConfigOption::ProcessHADR()
+{
+    // Process HADR option
+    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");
+    } 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.);
+    }
+}
+
+
+
+void TFlukaConfigOption::ProcessMULS()
+{
+    // Process MULS option
+    fprintf(fgFile,"*\n* --- MULS --- Muliple Scattering. Flag = %5d\n", fProcessFlag[kMULS]);
+    //
+    // Multiple scattering
+    //
+    if (fProcessFlag[kMULS] > 0) {
+       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);
+    }
+}
+
+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]);
+    //
+    // Ionisation energy loss
+    //
+    //
+    // Impose consistency
+    
+    if (fProcessFlag[kLOSS] == 1 || fProcessFlag[kLOSS] == 3) {
+       fProcessFlag[kDRAY] = 1;
+    } else if (fProcessFlag[kLOSS] == 2) {
+       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 (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;
+    Float_t stra = (Float_t) fProcessFlag[kSTRA];
+    
+
+    if (fProcessFlag[kLOSS] == 1 || fProcessFlag[kLOSS] == 3) {
+//
+// 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.);       
+    } 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.);      
+    }  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.);      
+    }
+}
+
+
+void TFlukaConfigOption::ProcessCUTGAM()
+{
+// Cut on gammas
+//
+    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.);
+    } 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.);
+       }
+    }
+}
+
+void TFlukaConfigOption::ProcessCUTELE()
+{
+// Cut on e+/e-
+//
+    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.);
+    } 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.);
+       }
+    }
+}
+
+void TFlukaConfigOption::ProcessCUTNEU()
+{
+    // Cut on neutral hadrons
+    fprintf(fgFile,"*\n*Cut for neutal hadrons. CUTNEU = %13.4g\n", fCutValue[kCUTNEU]);
+    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", -cut,  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);
+    } else {
+       printf("Cuts on neutral hadrons material by material not yet implemented !\n");
+    }
+}
+
+void TFlukaConfigOption::ProcessCUTHAD()
+{
+    // Cut on charged hadrons
+    fprintf(fgFile,"*\n*Cut for charge hadrons. CUTHAD = %13.4g\n", 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);
+    } else {
+       printf("Cuts on charged hadrons material by material not yet implemented !\n");
+    }
+}
+
+void TFlukaConfigOption::ProcessCUTMUO()
+{
+    // Cut on muons
+    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);
+    } else {
+       printf("Cuts on muons material by material not yet implemented !\n");
+    }
+    
+    
+}
+
+void TFlukaConfigOption::ProcessTOFMAX()
+{
+    // Cut time of flight
+    Float_t cut = fCutValue[kTOFMAX];
+    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);
 }
index 3153d58aa446e33b5e0f023f58e0249304543fb4..1b0fe05e3516bea3630a3a30026c57f7afec727e 100644 (file)
 //                                                                           //
 ///////////////////////////////////////////////////////////////////////////////
 //
-//
-//
 
 
-#include <TNamed.h>
+//
 
+#include <TObject.h>
 
-class TFlukaConfigOption : public TNamed
+typedef enum {kDCAY, kPAIR, kCOMP, kPHOT, kPFIS, kDRAY, kANNI, kBREM,
+             kMUNU, kCKOV, kHADR, kLOSS, kMULS, kRAYL, kSTRA} FlukaProcessOption_t;
+typedef enum {kCUTGAM, kCUTELE, kCUTNEU, kCUTHAD, kCUTMUO, kBCUTE, kBCUTM, kDCUTE, kDCUTM, kPPCUTM, kTOFMAX}  FlukaCutOption_t;
+class TFlukaMCGeometry;
+
+class TFlukaConfigOption : public TObject
 {
 public:
     // Constructors
     TFlukaConfigOption();
-    TFlukaConfigOption(const char* cutName, Double_t cut);
-    TFlukaConfigOption(const char* cutName, Double_t cut, Int_t imed);
-    TFlukaConfigOption(const char* procName, Int_t flag);
-    TFlukaConfigOption(const char* procName, Int_t flag, Int_t imed);
+    TFlukaConfigOption(Int_t imed);
     // Getters
-    Double_t Cut()  const                {return fCutValue;}
-    Int_t    Flag() const                {return fProcessFlag;}
-    Int_t    Medium()                    {return fMedium;}    
-    Bool_t   HasMediumAssigned()         {return  (fMedium > -1);}
+    Double_t Cut(FlukaCutOption_t i)      const {return fCutValue[i];}
+    Int_t    Flag(FlukaProcessOption_t i) const {return fProcessFlag[i];}
+    Int_t    Medium()                     const {return fMedium;}    
     // Setters
-    void     SetCut(Double_t val)        {fCutValue     =   val;}
-    void     SetFlag(Int_t val)          {fProcessFlag  =   val;}
-    void     SetMedium(Int_t imed)       {fMedium       =   imed;}    
-
+    void     SetCut(const char* flagname, Double_t val);
+    void     SetProcess(const char* flagname, Int_t flagValue);
+    void     SetMedium(Int_t imed)       {fMedium = imed;}
+    //
+    void     WriteFlukaInputCards();
+    void     ProcessDCAY();
+    void     ProcessPAIR();
+    void     ProcessBREM();
+    void     ProcessCOMP();
+    void     ProcessPHOT();
+    void     ProcessANNI();
+    void     ProcessPFIS();
+    void     ProcessMUNU();
+    void     ProcessRAYL();
+    void     ProcessCKOV();
+    void     ProcessHADR();
+    void     ProcessMULS();
+    void     ProcessLOSS();
+    void     ProcessSTRA();
+    
+    
+    void     ProcessCUTGAM();
+    void     ProcessCUTELE();
+    void     ProcessCUTNEU();
+    void     ProcessCUTHAD();
+    void     ProcessCUTMUO();
+    void     ProcessTOFMAX();
+    
+    //
+    static void SetStaticInfo(FILE* file, Float_t matMin, Float_t matMax, TFlukaMCGeometry* geom)
+       {fgFile = file; fgMatMin = matMin; fgMatMax = matMax; fgGeom = geom;}
  protected:
-    Double_t fCutValue;                // User cut
-    Int_t    fProcessFlag;             // User flag assigned to processes
+    Double_t fCutValue[11];            // User cut
+    Int_t    fProcessFlag[15];         // User flag assigned to processes
     Int_t    fMedium;                  // Materials assigned to user settings
+    Float_t  fCMatMin;                 // Minimum material number used for current card 
+    Float_t  fCMatMax;                 // Maximum material number used for current card
+    
+    // static
+    static Double_t  fgDCutValue[11];     // User default cut
+    static Int_t     fgDProcessFlag[15];  // User default flag assigned to processes
+    static Float_t   fgMatMin;            // Minimum material number 
+    static Float_t   fgMatMax;            // Maximum meterial number
+    static FILE*     fgFile;              // Output file
+    static TFlukaMCGeometry* fgGeom;      // Pointer to geometry     
     ClassDef(TFlukaConfigOption, 1)    // Fluka Configuration Option
 };