- Creation of peg files taking into account cuts on photon and electron energies
authormorsch <morsch@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 6 Jan 2005 07:17:59 +0000 (07:17 +0000)
committermorsch <morsch@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 6 Jan 2005 07:17:59 +0000 (07:17 +0000)
- Creation of pemf file separated from writing of material definitions.

TFluka/TFluka.cxx
TFluka/TFlukaMCGeometry.cxx
TFluka/TFlukaMCGeometry.h

index c6c722a..5772068 100644 (file)
@@ -168,8 +168,11 @@ TFluka::~TFluka() {
        fUserConfig->Delete();
        delete fUserConfig;
     }
-
-
+    
+    if (fUserScore) {
+       fUserScore->Delete();
+       delete fUserScore;
+    }
 }
 
 //
@@ -188,14 +191,16 @@ void TFluka::Init() {
     gGeoManager->SetTopVolume(top);
     gGeoManager->CloseGeometry("di");
     gGeoManager->DefaultColors();  // to be removed
+    
+    // Now we have TGeo geometry created and we have to patch FlukaVmc.inp
+    // with the material mapping file FlukaMat.inp
     fNVolumes = fGeom->NofVolumes();
     fGeom->CreateFlukaMatFile("flukaMat.inp");   
     if (fVerbosityLevel >=3) {
        printf("== Number of volumes: %i\n ==", fNVolumes);
        cout << "\t* InitPhysics() - Prepare input file to be called" << endl; 
-    }   
-    // now we have TGeo geometry created and we have to patch FlukaVmc.inp
-    // with the material mapping file FlukaMat.inp
+    }
 }
 
 
@@ -219,8 +224,34 @@ void TFluka::BuildPhysics() {
     
     if (fVerbosityLevel >=3)
        cout << "==> TFluka::BuildPhysics() called." << endl;
-// Prepare input file with the current physics settings
+
+    
+    if (fVerbosityLevel >=3) {
+       TList *medlist = gGeoManager->GetListOfMedia();
+       TIter next(medlist);
+       TGeoMedium*   med = 0x0;
+       TGeoMaterial* mat = 0x0;
+       Int_t ic = 0;
+       
+       while((med = (TGeoMedium*)next()))
+       {
+           mat = med->GetMaterial();
+           printf("Medium %5d %12s %5d %5d\n", ic, (med->GetName()), med->GetId(), mat->GetIndex());
+           ic++;
+       }
+    }
+    
+    //
+    // At this stage we have the information on materials and cuts available.
+    // Now create the pemf file
+    
+    if (fGeneratePemf) fGeom->CreatePemfFile();
+    
+    //
+    // Prepare input file with the current physics settings
+    
     InitPhysics(); 
+    
     cout << "\t* InitPhysics() - Prepare input file was called" << endl; 
     
     if (fVerbosityLevel >=2)
@@ -231,6 +262,7 @@ void TFluka::BuildPhysics() {
     if (fVerbosityLevel >=2)
        cout << "\t* Opening file " << fInputFileName << endl;
     const char* fname = fInputFileName;
+    
     fluka_openinp(lunin, PASSCHARA(fname));
     
     if (fVerbosityLevel >=2)
@@ -246,7 +278,6 @@ void TFluka::BuildPhysics() {
     if (fVerbosityLevel >=3)
        cout << "<== TFluka::Init() called." << endl;
     
-    
     if (fVerbosityLevel >=3)
        cout << "<== TFluka::BuildPhysics() called." << endl;
 }  
@@ -558,7 +589,7 @@ void TFluka::Gstpar(Int_t itmed, const char* param, Double_t parval) {
 //
 //
 // Check if material is used    
-   if (fVerbosityLevel >=3) 
+   if (fVerbosityLevel >= 3) 
        printf("Gstpar called with %6d %5s %12.4e %6d\n", itmed, param, parval, fGeom->GetFlukaMaterial(itmed));
    Int_t* reglist;
    Int_t nreg;
index 2f01abc..168df78 100644 (file)
@@ -27,6 +27,7 @@
 #include "TCallf77.h"
 #include "TFluka.h"
 #include "TFlukaMCGeometry.h"
+#include "TFlukaConfigOption.h"
 #include "TGeoManager.h" 
 #include "TGeoVolume.h" 
 #include "TObjString.h"
@@ -617,7 +618,8 @@ void TFlukaMCGeometry::CreateFlukaMatFile(const char *fname)
       fMatNames->Add(objstr);
       nfmater++;
    }
-   Int_t indmat = nfmater;
+   
+   fIndmat = nfmater;
 //   TGeoMedium *med;
    // Adjust material names and add them to FLUKA list
    for (i=0; i<nmater; i++) {
@@ -801,56 +803,65 @@ void TFlukaMCGeometry::CreateFlukaMatFile(const char *fname)
        if (gSystem->AccessPathName("FlukaVmc.pemf")) Fatal("CreateFlukaMatFile", "No pemf file in working directory");
        return;
    }
-   
-   // Write peg files
-   char number[20];
-   Int_t countMatOK = 0;
-   Int_t countElemError = 0;
-   Int_t countNoStern = 0;
-   Int_t countMixError = 0;
-   Int_t countGas = 0;
- //  Int_t countGasError = 0;
-   Int_t countPemfError = 0;
-   for (i=indmat; i<nfmater; i++) {
-      mat = (TGeoMaterial*)fMatList->At(i);
-      if (!mat->IsUsed()) continue;
-      sname = "mat";
-      sprintf(number, "%d", i);
-      sname.Append(number);
-      cout << endl;
-      cout << endl;
-      cout << "******************************************************************************" << endl;
-      cout << "******************************************************************************" << endl;
-      cout << endl;
-      WritePegFile(i, &countNoStern, &countElemError, &countMixError, &countGas);
-      sname.Prepend("$FLUPRO/pemf/rpemf peg/");
-      gSystem->Exec(sname.Data());
-      
-      // check if the pemf file was created
-      TString sname = Form("peg/mat%d.pemf", i);
-      ifstream in( sname.Data() );
-      if ( in ) {
-         countMatOK++;
-         in.close();
-      }
-      else {
-         cout << "ERROR Fail to create the pemf file " << sname << endl;
-         countPemfError++; 
-      }
-   }
-   cout << "Materials (pemf created)   " << countMatOK         << endl;
-   cout << "Not Sternheimer par. found  " << countNoStern   << endl;
-   cout << "Elements with error definitions (Z not integer)  " << countElemError      << endl;
-   cout << "Mixtures with error definitions (Z not integer) " << countMixError  << endl;
-   cout << "Posible Gas (rho < 0.01) " << countGas           << endl;
-  // cout << "Posible Gas (without pressure information) " << countGasError           << endl;
+}
+
+void TFlukaMCGeometry::CreatePemfFile()
+{
+    //
+    // Steering routine to write and process peg files producing the pemf input 
+    //
+    char number[20];
+    Int_t countMatOK     = 0;
+    Int_t countElemError = 0;
+    Int_t countNoStern   = 0;
+    Int_t countMixError  = 0;
+    Int_t countGas       = 0;
+    Int_t countPemfError = 0;
+    Int_t i;
+    TGeoMaterial* mat = 0x0;
+    TString sname;
+    
+    for (i = fIndmat; i < fLastMaterial - 2; i++) {
+       printf("Write Peg Files %d\n", i);
+       
+       mat = (TGeoMaterial*)fMatList->At(i);
+       if (!mat->IsUsed()) continue;
+       sname = "mat";
+       sprintf(number, "%d", i);
+       sname.Append(number);
+       cout << endl;
+       cout << endl;
+       cout << "******************************************************************************" << endl;
+       cout << "******************************************************************************" << endl;
+       cout << endl;
+       WritePegFile(i, &countNoStern, &countElemError, &countMixError, &countGas);
+       sname.Prepend("$FLUPRO/pemf/rpemf peg/");
+       gSystem->Exec(sname.Data());
+       
+       // check if the pemf file was created
+       TString sname = Form("peg/mat%d.pemf", i);
+       ifstream in( sname.Data() );
+       if ( in ) {
+           countMatOK++;
+           in.close();
+       } else {
+           cout << "ERROR Fail to create the pemf file " << sname << endl;
+           countPemfError++; 
+       }
+    }
+    cout << "Materials (pemf created)   " << countMatOK         << endl;
+    cout << "Not Sternheimer par. found  " << countNoStern   << endl;
+    cout << "Elements with error definitions (Z not integer)  " << countElemError      << endl;
+    cout << "Mixtures with error definitions (Z not integer) " << countMixError  << endl;
+    cout << "Posible Gas (rho < 0.01) " << countGas           << endl;
+    // cout << "Posible Gas (without pressure information) " << countGasError           << endl;
     cout << "Pemf files Error    " << countPemfError     << endl;
-   cout << endl << endl;
-   
-   sname = "cat peg/*.pemf > peg/FlukaVmc.pemf";         
-   gSystem->Exec(sname.Data());
-   sname = "mv peg/FlukaVmc.pemf FlukaVmc.pemf";
-   gSystem->Exec(sname.Data());
+    cout << endl << endl;
+    
+    sname = "cat peg/*.pemf > peg/FlukaVmc.pemf";         
+    gSystem->Exec(sname.Data());
+    sname = "mv peg/FlukaVmc.pemf FlukaVmc.pemf";
+    gSystem->Exec(sname.Data());
 }
 
 //_____________________________________________________________________________
@@ -996,8 +1007,35 @@ void TFlukaMCGeometry::WritePegFile(Int_t imat, Int_t *NoStern, Int_t *ElemError
       delete [] wt;
    }
    
+   Double_t ue = 3000000.; // [MeV]
+   Double_t up = 3000000.; // [MeV]
+   Double_t ae = -1.;
+   Double_t ap = -1.;
+   
+   
+   TObjArray* cutList = ((TFluka*) gMC)->GetListOfUserConfigs();
+   TIter next(cutList);
+   TFlukaConfigOption* proc;
+   
+   while((proc = (TFlukaConfigOption*)next()))
+   { 
+       if (proc->Medium() == mat->GetIndex()) {
+          ap = proc->Cut(kCUTGAM);
+          ae = proc->Cut(kCUTELE);
+          if (ap == -1.) ap =  TFlukaConfigOption::DefaultCut(kCUTGAM);
+          if (ae == -1.) ae =  TFlukaConfigOption::DefaultCut(kCUTELE);
+          break;
+       }
+   }
+
+   if (ap == -1.) ap =  TFlukaConfigOption::DefaultCut(kCUTGAM);
+   if (ae == -1.) ae =  TFlukaConfigOption::DefaultCut(kCUTELE);
+
+   ap *= 1000.;                         // [MeV]
+   ae  = (ae + 0.00051099906) * 1000.;  // [MeV]
+   
    out << "ENER" << endl;
-   out << " $INP AE=0.56099906, UE=3000000., AP=.03, UP=3000000. $END" << endl;
+   out << " $INP AE=" << ae << ", UE=" << ue <<", AP=" << ap << ", UP=" << up << " $END" << endl;
    out << "PWLF" << endl;
    out << " $INP NALE=300, NALG=400, NALR=100 $END" << endl;
    out << "DECK" << endl;
@@ -1023,7 +1061,7 @@ Double_t * TFlukaMCGeometry::GetISSB(Double_t rho, Int_t nElem, Double_t *zelem,
       Double_t    wt[20];             //[nelems] weight fraction
       Double_t    density;            // g/cm3
       Double_t    iev;                // Average Ion potential (eV)
-                                   // ****   Sternheimer parameters  ****
+                                      // ****   Sternheimer parameters  ****
       Double_t    cbar;               // CBAR
       Double_t    x0;                 // X0
       Double_t    x1;                 // X1
index 88ef81d..d23a32f 100644 (file)
@@ -37,9 +37,9 @@ class TFlukaMCGeometry :public TNamed {
     Int_t         GetFlukaMaterial(Int_t imed) const;
     Int_t         GetLastMaterialIndex() const {return fLastMaterial;}
     virtual Int_t NofVolumes() const;
-
    // FLUKA specific methods
     void          CreateFlukaMatFile(const char *fname=0);
+    void          CreatePemfFile();
     void          PrintHeader(ofstream &out, const char *text) const;
     Bool_t        IsDebugging() const {return fDebug;}
     void          SetDebugMode(Bool_t flag=kTRUE) {fDebug = flag;}
@@ -74,6 +74,7 @@ class TFlukaMCGeometry :public TNamed {
     Int_t        fNextRegion;             // next region number
     Int_t        fNextLattice;            // next lattice history
     Int_t       *fRegionList;             //! region list matching a given medium number
+    Int_t        fIndmat;                 // material index where pemf file creation starts
     TObjArray   *fMatList;                //! material list as known by FLUKA
     TObjArray   *fMatNames;               //! list of FLUKA material names
   ClassDef(TFlukaMCGeometry,1)  //Virtual MonteCarlo Interface