Pemf generation option is obsolete. Removed.
authormorsch <morsch@f7af4fe6-9843-0410-8265-dc069ae4e863>
Mon, 15 Oct 2007 07:48:14 +0000 (07:48 +0000)
committermorsch <morsch@f7af4fe6-9843-0410-8265-dc069ae4e863>
Mon, 15 Oct 2007 07:48:14 +0000 (07:48 +0000)
TFluka/TFluka.cxx
TFluka/TFluka.h
TFluka/TFlukaMCGeometry.cxx
TFluka/TFlukaMCGeometry.h

index 1c3a90f..80faf52 100644 (file)
@@ -133,7 +133,6 @@ TFluka::TFluka()
    fTrackIsExiting(kFALSE),
    fTrackIsNew(kFALSE),
    fFieldFlag(kTRUE),
-   fGeneratePemf(kFALSE),
    fDummyBoundary(kFALSE),
    fStopped(kFALSE),
    fStopEvent(kFALSE),
@@ -171,7 +170,6 @@ TFluka::TFluka(const char *title, Int_t verbosity, Bool_t isRootGeometrySupporte
    fTrackIsExiting(kFALSE),
    fTrackIsNew(kFALSE),
    fFieldFlag(kTRUE),
-   fGeneratePemf(kFALSE),
    fDummyBoundary(kFALSE),
    fStopped(kFALSE),
    fStopEvent(kFALSE),
@@ -297,13 +295,7 @@ void TFluka::BuildPhysics() {
         }
     }
     
-    //
-    // 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(); 
index 3ef6f42..7909d0c 100644 (file)
@@ -378,10 +378,6 @@ class TFluka : public TVirtualMC {
   Int_t  GetDummyBoundary() const {return fDummyBoundary;}
   Bool_t IsDummyBoundary() const {return (fDummyBoundary==0)?kFALSE:kTRUE;}
   
-
-  void   SetGeneratePemf(Bool_t flag = kTRUE) {fGeneratePemf = flag;}
-  Bool_t IsGeneratePemf() const {return fGeneratePemf;}
-  
   void   EnableField(Bool_t flag=kTRUE) {fFieldFlag = flag;}
   Bool_t IsFieldEnabled() const {return fFieldFlag;}
   
@@ -428,7 +424,6 @@ class TFluka : public TVirtualMC {
   Bool_t   fTrackIsExiting;       // Flag for track exiting  
   Bool_t   fTrackIsNew;           // Flag for new track
   Bool_t   fFieldFlag;            // Flag for magnetic field
-  Bool_t   fGeneratePemf;         // Flag for automatic .pemf generation
   Int_t    fDummyBoundary;        // Flag for crossing dummy boundaries
   Bool_t   fStopped;              // Flag for stopping 
   Bool_t   fStopEvent;            // Flag for stopped event
index 1199d3c..295717c 100644 (file)
@@ -810,387 +810,6 @@ void TFlukaMCGeometry::CreateFlukaMatFile(const char *fname)
    fLastMaterial = nfmater+2;
 }
 
-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());
-}
-
-//_____________________________________________________________________________
-void TFlukaMCGeometry::WritePegFile(Int_t imat, Int_t *NoStern, Int_t *ElemError,
-                       Int_t *MixError, Int_t *countGas) const
-{
-   // Write the .peg file for one material
-   
-   TGeoMaterial *mat = (TGeoMaterial*)fMatList->At(imat);
-   TString name = ((TObjString*)fMatNames->At(imat))->GetString();
-   TString line;
-   char number[20];
-   TGeoElement *elem = mat->GetElement();
-   name = name.Strip();
-   TString sname = "mat";
-   sprintf(number, "%d", imat);
-   sname.Append(number);
-   sname.Append(".peg");
-   sname.Prepend("peg/");
-   ofstream out;
-   out.open(sname.Data(), ios::out);
-   if (!out.good()) return;
-   Double_t dens = mat->GetDensity();
-   TGeoMixture *mix = 0;
-   Int_t nel = 1;
-   Int_t i;
-   if (mat->IsMixture()) {
-      mix = (TGeoMixture*)mat;
-      nel = mix->GetNelements();
-   } 
-     
-   if (nel==1) {
-      cout  << "( Element ) " << name << "  Z=" << mat->GetZ() << " Rho " << mat->GetDensity() << endl;
-
-      Double_t zel = mat->GetZ();
-      if( (zel-Int_t(zel))>0.001 || zel < 1 ) {
-         cout << " ERROR: A Element with not integer Z=" << zel << endl;
-         cout << endl;
-         (*ElemError)++;
-         return;
-      }
-      
-      out << "ELEM" << endl;
-      out << " &INP IRAYL=1, RHO=" << dens << ", " << endl;
-      
-      // check for the Sternheimer parameters
-      Double_t *issb_parm = GetISSB( mat->GetDensity(), 1, &zel, 0 );
-      if( issb_parm[0] > 0 && issb_parm[1] > 0 ) {
-         cout << "Sternheimer parameters found" << endl;
-         out << ", ISSB=1, IEV=" << issb_parm[0] << ", CBAR=" << issb_parm[1]
-             << ", X0=" << issb_parm[2] << "," << endl;
-         out << "X1=" <<issb_parm[3] <<", AFACT="<<issb_parm[4] <<", SK="
-             << issb_parm[5] << ", DELTA0=" << issb_parm[6];
-      } 
-      else {
-         cout << "WARNING: Strange element, Sternheimer parameters  not found" << endl;
-        (*NoStern)++;
-      }
-
-      if (dens<0.01) {
-        (*countGas)++;
-        out << " GASP=1." << endl;
-      }
-      
-      out << " &END" <<  endl;
-      out << name.Data() << endl;
-      out << elem->GetName() << endl;
-      
-   } 
-   else {
-   
-      cout  << "( Mixture )  " << name << "  Rho " << dens << " nElem " << nel << endl;
-    
-      Double_t *zt = new Double_t[nel];
-      Double_t *wt = new Double_t[nel];
-      for (int j=0; j<nel; j++) {
-         zt[j] = (mix->GetZmixt())[j];
-         wt[j] = (mix->GetWmixt())[j];
-         if( (zt[j]-Int_t(zt[j])) > 0.001 || zt[j] < 1 ) {
-            cout << "ERROR Mixture " << name << " with an element with not integer Z=" << zt[j] << endl;
-            cout << endl;
-            (*MixError)++;
-            // just continue since the mixtures are not patch, 
-            // but the final release should include the return   
-            //  return;         
-         }
-      }
-      Double_t *issb_parm = GetISSB( mat->GetDensity(), nel, zt, wt );
-      out << "MIXT" << endl;
-      out << " &INP IRAYL=1, NE=" << nel << ", RHOZ=" << wt[0] << ",";
-      line = Form(" &INP IRAYL=1, NE=%d RHOZ=%g", nel, wt[0]);
-      for(int j=1; j<nel; j++) {
-         out << " " << wt[j] << ",";
-         line += Form(" %g,", wt[j] );
-         if( line.Length() > 60 ) { out << endl; line = ""; }
-      }
-      out << " RHO=" << mat->GetDensity() << ", ";
-      line += Form(" RHO=%g, ", mat->GetDensity());
-      if( line.Length() > 60 ) { out << endl; line = ""; }
-      
-      if( issb_parm[0] > 0 && issb_parm[1] > 0 ) {
-         cout << "Sternheimer parameters found" << endl;
-         out << " ISSB=1, IEV=" << issb_parm[0] << ",";
-         line += Form(" ISSB=1, IEV=%g,", issb_parm[0]);
-         if( line.Length() > 60 ) { out << endl; line = "";  }
-         out << " CBAR=" << issb_parm[1] << ",";
-         line += Form(" CBAR=%g,",issb_parm[1]);
-         if( line.Length() > 60 ) { out << endl; line = "";  }
-         out << " X0=" << issb_parm[2] << ",";
-         line += Form(" X0=%g,", issb_parm[2]);
-         if( line.Length() > 60 ) { out << endl; line = "";  }
-         out << " X1=" << issb_parm[3] << ",";
-         line += Form(" X1=%g,", issb_parm[3]);
-         if( line.Length() > 60 ) { out << endl; line = "";  }
-         out << " AFACT="<< issb_parm[4] << ",";
-         line += Form(" AFACT=%g,", issb_parm[4]);
-         if( line.Length() > 60 ) { out << endl; line = "";  }
-         out << " SK=" << issb_parm[5] << ",";
-         line += Form(" SK=%g,", issb_parm[5]);
-         if( line.Length() > 60 ) { out << endl; line = "";  }
-      }
-      else {
-         cout << "Sternheimer parameters  not found" << endl;
-         (*NoStern)++;
-      }
-      
-      if (dens<0.01){
-         (*countGas)++;
-         out << " GASP=1." << endl;
-      }
-      
-      out << " &END" <<  endl;
-      out << name.Data() << endl;
-      for (i=0; i<nel; i++) {
-         elem = mix->GetElement(i);
-         line = elem->GetName();
-         if (line.Length()==1) line.Append(" ");
-         out << line.Data() << " ";
-      }
-      out << endl;
-      
-      delete [] zt;
-      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=" << 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;
-   out << " $INP $END" << endl;
-   out << "TEST" << endl;
-   out << " $INP $END" << endl;
-   out.close();
-}
-
-Double_t * TFlukaMCGeometry::GetISSB(Double_t rho, Int_t nElem, Double_t *zelem, Double_t *welem ) const
-{
-   // Read the density effect parameters
-   // from R.M. Sternheimer et al. Atomic Data
-   // and Nuclear Data Tables, Vol. 30 No. 2
-   //
-   // return the parameters if the element/mixture match with one of the list
-   // otherwise returns the parameters set to 0
-   
-   struct sternheimerData {
-      TString     longname;           // element/mixture name
-      Int_t       nelems;             // number of constituents N
-      Int_t       Z[20];              //[nelems] Z
-      Double_t    wt[20];             //[nelems] weight fraction
-      Double_t    density;            // g/cm3
-      Double_t    iev;                // Average Ion potential (eV)
-                                      // ****   Sternheimer parameters  ****
-      Double_t    cbar;               // CBAR
-      Double_t    x0;                 // X0
-      Double_t    x1;                 // X1
-      Double_t    afact;              // AFACT
-      Double_t    sk;                 // SK
-      Double_t    delta0;             // DELTA0
-
-      sternheimerData():
-         longname(""), nelems(0), density(0), iev(0), cbar(0),
-         x0(0), x1(0), afact(0), sk(0), delta0(0) {}
-   };
-   
-   TString     shortname;
-   TString     formula;
-   Int_t       num;
-   char        state;
-   
-   static Double_t parameters[7];
-   memset( parameters, 0, sizeof(Double_t) );
-
-   static sternheimerData sternDataArray[300];
-   static Bool_t isFileRead = kFALSE;
-   
-   // Read the data file if is needed
-   if( isFileRead == kFALSE ) {
-      TString sSternheimerInp = getenv("ALICE_ROOT");
-      sSternheimerInp +="/TFluka/input/Sternheimer.data";
-   
-      ifstream in(sSternheimerInp);
-      char line[100];
-      in.getline(line, 100);   
-      in.getline(line, 100);   
-      in.getline(line, 100);   
-      in.getline(line, 100);   
-      in.getline(line, 100);   
-      in.getline(line, 100);   
-      
-      
-      Int_t is = 0;
-      while( !in.eof() ) {
-         in >> shortname >> num     >> sternDataArray[is].nelems 
-            >> sternDataArray[is].longname  >> formula >> state;
-         if( in.eof() ) break;
-         for(int i=0; i<sternDataArray[is].nelems; i++) {
-            in >> sternDataArray[is].Z[i] >> sternDataArray[is].wt[i]; 
-         }
-         in >> sternDataArray[is].density; 
-         in >> sternDataArray[is].iev; 
-         in >> sternDataArray[is].cbar; 
-         in >> sternDataArray[is].x0; 
-         in >> sternDataArray[is].x1; 
-         in >> sternDataArray[is].afact; 
-         in >> sternDataArray[is].sk;
-         if( sternDataArray[is].nelems == 1 ) in >> sternDataArray[is].delta0;
-         is++;
-      }
-      isFileRead = kTRUE;
-      in.close();
-   }   
-   
-   Int_t is = 0;
-   while( is < 280 ) {
-   
-      // check for elements
-      if( sternDataArray[is].nelems == 1 && nElem == 1
-          && sternDataArray[is].Z[0] == Int_t(*zelem)
-          && TMath::Abs( (sternDataArray[is].density - rho)/sternDataArray[is].density ) < 0.1 ) {
-         cout << sternDataArray[is].longname << "   #elems:" <<  sternDataArray[is].nelems << "  Rho:" 
-              << sternDataArray[is].density << endl;
-         cout << sternDataArray[is].iev   << " " 
-              << sternDataArray[is].cbar  << " " 
-              << sternDataArray[is].x0    << " " 
-              << sternDataArray[is].x1    << " " 
-              << sternDataArray[is].afact << " " 
-              << sternDataArray[is].sk    << " " 
-              << sternDataArray[is].delta0 << endl;
-         
-         parameters[0] = sternDataArray[is].iev;
-         parameters[1] = sternDataArray[is].cbar;
-         parameters[2] = sternDataArray[is].x0;
-         parameters[3] = sternDataArray[is].x1;
-         parameters[4] = sternDataArray[is].afact;
-         parameters[5] = sternDataArray[is].sk;
-         parameters[6] = sternDataArray[is].delta0;
-         return parameters;        
-      }
-      
-      // check for mixture
-      int nmatch = 0;
-      if( sternDataArray[is].nelems > 1 && sternDataArray[is].nelems == nElem ) {
-         for(int j=0; j<sternDataArray[is].nelems; j++) {
-            if( sternDataArray[is].Z[j] == Int_t(zelem[j]) && 
-               TMath::Abs( (sternDataArray[is].wt[j] - welem[j])/sternDataArray[is].wt[j] ) < 0.1 )
-            nmatch++;            
-         }
-      }
-
-      if( sternDataArray[is].nelems > 1 && 
-          TMath::Abs( (sternDataArray[is].density - rho)/sternDataArray[is].density ) < 0.1 
-          && nmatch == sternDataArray[is].nelems ) {
-         cout << sternDataArray[is].longname << "   #elem:" <<  sternDataArray[is].nelems << "  Rho:" 
-              << sternDataArray[is].density << endl;
-         cout << sternDataArray[is].iev   << " " 
-              << sternDataArray[is].cbar  << " " 
-              << sternDataArray[is].x0    << " " 
-              << sternDataArray[is].x1    << " " 
-              << sternDataArray[is].afact << " " 
-              << sternDataArray[is].sk    << " " 
-              << sternDataArray[is].delta0 << endl;
-
-         parameters[0] = sternDataArray[is].iev;
-         parameters[1] = sternDataArray[is].cbar;
-         parameters[2] = sternDataArray[is].x0;
-         parameters[3] = sternDataArray[is].x1;
-         parameters[4] = sternDataArray[is].afact;
-         parameters[5] = sternDataArray[is].sk;
-         parameters[6] = 0;
-         return parameters;        
-      }
-      is++; 
-   }   
-   return parameters;        
-}
-
 //_____________________________________________________________________________
 void TFlukaMCGeometry::PrintHeader(ofstream &out, const char *text) const
 {
index d7dac6a..768530a 100644 (file)
@@ -40,7 +40,6 @@ class TFlukaMCGeometry :public TNamed {
     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;}
@@ -60,12 +59,8 @@ class TFlukaMCGeometry :public TNamed {
     void          ToFlukaString(TString &str) const;
     void          FlukaMatName(TString &str) const;
     Int_t         GetPredefinedMaterialId(Int_t z) const;
-    void          WritePegFile(Int_t imat, Int_t *NoStern, Int_t *ElemError,
-                       Int_t *MixError, Int_t *countGas) const;
-    Double_t *    GetISSB(Double_t rho, Int_t nElem, Double_t *zelem, Double_t *welem ) const;
-
-    Double_t* CreateDoubleArray(Float_t* array, Int_t size) const;
-    void     Vname(const char *name, char *vname) const;
+    Double_t*     CreateDoubleArray(Float_t* array, Int_t size) const;
+    void          Vname(const char *name, char *vname) const;
    
   private:
     // Copy constructor and operator= declared but not implemented (-Weff++ flag)