- Correct setting of FUDGEM parameter.

[u/mrichter/AliRoot.git] / TFluka / TFlukaConfigOption.cxx

diff --git a/TFluka/TFlukaConfigOption.cxx b/TFluka/TFlukaConfigOption.cxx
index e8d7bb53ebfc3598b5394dcaee5ac219713091e0..d6dde3fc6dfc07e0b0389c6fd3fdf79ee0373fab 100644 (file)
--- a/TFluka/TFlukaConfigOption.cxx
+++ b/TFluka/TFlukaConfigOption.cxx
@@ -21,9 +21,13 @@
 #include "TFlukaCerenkov.h"
 
 #include <TString.h>
 #include "TFlukaCerenkov.h"
 
 #include <TString.h>

+#include <TList.h>

 #include <TObjArray.h>
 #include <TVirtualMC.h>
 #include <TGeoMaterial.h>
 #include <TObjArray.h>
 #include <TVirtualMC.h>
 #include <TGeoMaterial.h>

+#include <TGeoMedium.h>
+#include <TGeoManager.h>
+#include <TGeoMedium.h>

 
 Float_t            TFlukaConfigOption::fgMatMin(-1.);
 Float_t            TFlukaConfigOption::fgMatMax(-1.);
 
 Float_t            TFlukaConfigOption::fgMatMin(-1.);
 Float_t            TFlukaConfigOption::fgMatMax(-1.);


@@ -97,6 +101,14 @@ void TFlukaConfigOption::WriteFlukaInputCards()
     // Write the FLUKA input cards for the set of process flags and cuts
     //
     //
     // Write the FLUKA input cards for the set of process flags and cuts
     //
     //

+    //
+    // Check if global option or medium specific
+
+    Bool_t mediumIsSensitive = kFALSE;
+    TGeoMedium*   med    = 0x0;
+    TGeoMedium*   medium = 0x0;    
+    TGeoMaterial* mat    = 0x0;
+

     if (fMedium != -1) {
        TFluka* fluka = (TFluka*) gMC;
        fMedium = fgGeom->GetFlukaMaterial(fMedium);
     if (fMedium != -1) {
        TFluka* fluka = (TFluka*) gMC;
        fMedium = fgGeom->GetFlukaMaterial(fMedium);


@@ -109,22 +121,35 @@ void TFlukaConfigOption::WriteFlukaInputCards()
            printf("Material not used !\n");
            return;
        }
            printf("Material not used !\n");
            return;
        }

-
+       //
+       // Find material

        TObjArray *matList = fluka->GetFlukaMaterials();
        Int_t nmaterial =  matList->GetEntriesFast();
        TObjArray *matList = fluka->GetFlukaMaterials();
        Int_t nmaterial =  matList->GetEntriesFast();

-       TGeoMaterial* material = 0;
+       fCMaterial = 0;

        for (Int_t im = 0; im < nmaterial; im++)
        {
        for (Int_t im = 0; im < nmaterial; im++)
        {

-           material = dynamic_cast<TGeoMaterial*> (matList->At(im));
-           Int_t idmat = material->GetIndex();
+           fCMaterial = dynamic_cast<TGeoMaterial*> (matList->At(im));
+           Int_t idmat = fCMaterial->GetIndex();

            if (idmat == fMedium) break;            
            if (idmat == fMedium) break;            

-       }
-       
-       
-//
-// Check if global option
+       } // materials
+        //
+       // Find medium
+       TList *medlist = gGeoManager->GetListOfMedia();
+       TIter next(medlist);
+       while((med = (TGeoMedium*)next()))
+       {
+           mat = med->GetMaterial();
+           if (mat->GetIndex() == fMedium) {
+               medium = med;
+               break;
+           }
+       } // media
+       //
+       // Check if sensitive
+       if (medium->GetParam(0) != 0.) mediumIsSensitive = kTRUE;
+

 
 

-       fprintf(fgFile,"*\n*Material specific process and cut settings for #%8d %s\n", fMedium, material->GetName());
+       fprintf(fgFile,"*\n*Material specific process and cut settings for #%8d %s\n", fMedium, fCMaterial->GetName());

        fCMatMin = fMedium;
        fCMatMax = fMedium;
     } else {
        fCMatMin = fMedium;
        fCMatMax = fMedium;
     } else {


@@ -135,7 +160,17 @@ void TFlukaConfigOption::WriteFlukaInputCards()
 
 //
 // Handle Process Flags 
 
 //
 // Handle Process Flags 

+//
+//
+//  First make sure that all cuts are taken into account
+    if (DefaultProcessFlag(kPAIR) > 0 && fProcessFlag[kPAIR] == -1 && (fCutValue[kCUTELE] >= 0. || fCutValue[kPPCUTM] >= 0.)) 
+       fProcessFlag[kPAIR] = DefaultProcessFlag(kPAIR);
+    if (DefaultProcessFlag(kBREM) > 0 && fProcessFlag[kBREM] == -1 && (fCutValue[kBCUTE]  >= 0. || fCutValue[kBCUTM] >= 0.)) 
+       fProcessFlag[kBREM] = DefaultProcessFlag(kBREM);
+    if (DefaultProcessFlag(kDRAY) > 0 && fProcessFlag[kDRAY] == -1 && (fCutValue[kDCUTE]  >= 0. || fCutValue[kDCUTM] >= 0.)) 
+       fProcessFlag[kDRAY] = DefaultProcessFlag(kDRAY);

 //    
 //    

+//

     if (fProcessFlag[kDCAY] != -1) ProcessDCAY();
     if (fProcessFlag[kPAIR] != -1) ProcessPAIR();
     if (fProcessFlag[kBREM] != -1) ProcessBREM();
     if (fProcessFlag[kDCAY] != -1) ProcessDCAY();
     if (fProcessFlag[kPAIR] != -1) ProcessPAIR();
     if (fProcessFlag[kBREM] != -1) ProcessBREM();


@@ -159,8 +194,13 @@ void TFlukaConfigOption::WriteFlukaInputCards()
     if (fCutValue[kCUTNEU] >= 0.) ProcessCUTNEU();
     if (fCutValue[kCUTHAD] >= 0.) ProcessCUTHAD();
     if (fCutValue[kCUTMUO] >= 0.) ProcessCUTMUO();
     if (fCutValue[kCUTNEU] >= 0.) ProcessCUTNEU();
     if (fCutValue[kCUTHAD] >= 0.) ProcessCUTHAD();
     if (fCutValue[kCUTMUO] >= 0.) ProcessCUTMUO();

-
+//
+//  Time of flight 

     if (fCutValue[kTOFMAX] >= 0.) ProcessTOFMAX();
     if (fCutValue[kTOFMAX] >= 0.) ProcessTOFMAX();

+
+//
+//  Tracking precission
+    if (mediumIsSensitive) ProcessSensitiveMedium();

 }
 
 void TFlukaConfigOption::ProcessDCAY()
 }
 
 void TFlukaConfigOption::ProcessDCAY()


@@ -543,8 +583,14 @@ void TFlukaConfigOption::ProcessCUTGAM()
            fprintf(fgFile,"EMFCUT    %10.4g%10.4g%10.1f%10.1f%10.1f%10.1f\n", 0.,fCutValue[kCUTGAM], 0., ireg, ireg, 1.);
        }
     }
            fprintf(fgFile,"EMFCUT    %10.4g%10.4g%10.1f%10.1f%10.1f%10.1f\n", 0.,fCutValue[kCUTGAM], 0., ireg, ireg, 1.);
        }
     }

+    
+    // Transport production cut used for pemf
+    //
+    //  FUDGEM paramter. The parameter takes into account th contribution of atomic electrons to multiple scattering.
+    //  For production and transport cut-offs larger than 100 keV it must be set = 1.0, while in the keV region it must be
+    Float_t parFudgem = (fCutValue[kCUTGAM] > 1.e-4)? 1.0 : 0.0 ;

     fprintf(fgFile,"EMFCUT    %10.4g%10.4g%10.1f%10.1f%10.1f%10.1fPROD-CUT\n", 
     fprintf(fgFile,"EMFCUT    %10.4g%10.4g%10.1f%10.1f%10.1f%10.1fPROD-CUT\n", 

-           0., fCutValue[kCUTGAM], 1., fCMatMin, fCMatMax, 1.);
+           0., fCutValue[kCUTGAM], parFudgem, fCMatMin, fCMatMax, 1.);

 }
 
 void TFlukaConfigOption::ProcessCUTELE()
 }
 
 void TFlukaConfigOption::ProcessCUTELE()


@@ -565,8 +611,13 @@ void TFlukaConfigOption::ProcessCUTELE()
            fprintf(fgFile,"EMFCUT    %10.4g%10.4g%10.1f%10.1f%10.1f%10.1f\n", -fCutValue[kCUTELE], 0., 0., ireg, ireg, 1.);
        }
     }
            fprintf(fgFile,"EMFCUT    %10.4g%10.4g%10.1f%10.1f%10.1f%10.1f\n", -fCutValue[kCUTELE], 0., 0., ireg, ireg, 1.);
        }
     }

+    // Transport production cut used for pemf
+    //
+    //  FUDGEM paramter. The parameter takes into account th contribution of atomic electrons to multiple scattering.
+    //  For production and transport cut-offs larger than 100 keV it must be set = 1.0, while in the keV region it must be
+    Float_t parFudgem = (fCutValue[kCUTELE] > 1.e-4)? 1.0 : 0.0;

     fprintf(fgFile,"EMFCUT    %10.4g%10.4g%10.1f%10.1f%10.1f%10.1fPROD-CUT\n", 
     fprintf(fgFile,"EMFCUT    %10.4g%10.4g%10.1f%10.1f%10.1f%10.1fPROD-CUT\n", 

-           -fCutValue[kCUTELE], 0., 1., fCMatMin, fCMatMax, 1.);
+           -fCutValue[kCUTELE], 0., parFudgem, fCMatMin, fCMatMax, 1.);

 }
 
 void TFlukaConfigOption::ProcessCUTNEU()
 }
 
 void TFlukaConfigOption::ProcessCUTNEU()


@@ -654,7 +705,9 @@ void TFlukaConfigOption::ProcessCUTNEU()
                  Float_t(groupCut), 73.0, 0.95, ireg, ireg, 1.);
        }
 
                  Float_t(groupCut), 73.0, 0.95, ireg, ireg, 1.);
        }
 

-       printf("Cuts on neutral hadrons material by material only implemented for low-energy neutrons !\n");
+       Warning("ProcessCUTNEU", 
+               "Material #%4d %s: Cut on neutral hadrons (Ekin > %9.3e) material by material only implemented for low-energy neutrons !\n", 
+               fMedium, fCMaterial->GetName(), cut);

     }
 }
 
     }
 }
 


@@ -662,8 +715,8 @@ void TFlukaConfigOption::ProcessCUTHAD()
 {
     // Cut on charged hadrons
     fprintf(fgFile,"*\n*Cut for charge hadrons. CUTHAD = %13.4g\n", fCutValue[kCUTHAD]);
 {
     // Cut on charged hadrons
     fprintf(fgFile,"*\n*Cut for charge hadrons. CUTHAD = %13.4g\n", fCutValue[kCUTHAD]);

+    Float_t cut = fCutValue[kCUTHAD];

     if (fMedium == -1) {
     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);
        // 1.0 = Proton
        // 2.0 = Antiproton
        fprintf(fgFile,"PART-THR  %10.4g%10.1f%10.1f\n", -cut,  1.0,  2.0);


@@ -694,7 +747,9 @@ void TFlukaConfigOption::ProcessCUTHAD()
        // 58.0 = AntiXi_c minus
        fprintf(fgFile,"PART-THR  %10.4g%10.1f%10.1f\n", -cut, 57.0, 58.0);
     } else {
        // 58.0 = AntiXi_c minus
        fprintf(fgFile,"PART-THR  %10.4g%10.1f%10.1f\n", -cut, 57.0, 58.0);
     } else {

-      printf("Cuts on charged hadrons material by material not yet implemented !\n"); 
+      Warning("ProcessCUTHAD", 
+             "Material #%4d %s: Cut on charged hadrons (Ekin > 9.3e) material by material not yet implemented !\n", 
+             fMedium, fCMaterial->GetName(), cut); 

     }
 }
 
     }
 }
 


@@ -706,7 +761,8 @@ void TFlukaConfigOption::ProcessCUTMUO()
     if (fMedium == -1) {
        fprintf(fgFile,"PART-THR  %10.4g%10.1f%10.1f\n",-cut, 10.0, 11.0);
     } else {
     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");
+       Warning("ProcessCUTMUO", "Material #%4d %s: Cut on muons (Ekin > %9.3e) material by material not yet implemented !\n", 
+               fMedium, fCMaterial->GetName(), cut);

     }
     
     
     }
     
     


@@ -719,3 +775,18 @@ void TFlukaConfigOption::ProcessTOFMAX()
     fprintf(fgFile,"*\n*Cut on time of flight. TOFMAX = %13.4g\n", fCutValue[kTOFMAX]);
     fprintf(fgFile,"TIME-CUT  %10.4g%10.1f%10.1f%10.1f%10.1f\n",cut*1.e9,0.,0.,-6.0,64.0);
 }
     fprintf(fgFile,"*\n*Cut on time of flight. TOFMAX = %13.4g\n", fCutValue[kTOFMAX]);
     fprintf(fgFile,"TIME-CUT  %10.4g%10.1f%10.1f%10.1f%10.1f\n",cut*1.e9,0.,0.,-6.0,64.0);
 }

+
+void  TFlukaConfigOption::ProcessSensitiveMedium()
+{
+    //
+    // Special options for sensitive media
+    //
+
+    fprintf(fgFile,"*\n*Options for sensitive medium \n");
+    //
+    // EMFFIX
+    fprintf(fgFile,"EMFFIX    %10.1f%10.3f%10.1f%10.1f%10.1f%10.1f\n", fCMatMin, 0.05, 0., 0., 0., 0.);
+    //
+    // FLUKAFIX
+    fprintf(fgFile,"FLUKAFIX  %10.3f                    %10.3f\n", 0.05, fCMatMin);
+}