[g4star] Miscellaneous updates to StarAgmlLib (#308)

* Add map of GEANT3 paramaters and cuts associated to volumes.

* AgMLExtension inherits from TGeoExtension, which allows user and/or
framework information to be attached to TGeoVolumes.  The VMC framework
utlizes the "framework" extensions, so we will utilize the "user"
extensions in AgML.

- Attach user-defined hit scoring functions to sensitive volumes.
  When the MC application's hit scoring is invoked in a given
  volume, the general hit scoring function can invoke any user
  functions stored on the volume.

- Store the number of volume branchings for each volume, allowing
  us to define a unique volume ID to each volume defined within
  an AgML module.

* Cleanup memory.

* General cleanup, and moving mixture component to protected.

* AgML: Introduce a namespace, and uncomment a few trig functions

AgMath provides math functions used by AgML codes. A few of these are
functions defined in FORtran. A few defines have been made to prevent
AgML variables from being mistaken for c-math routines.

* Attach AgMLExtensions to the modules

* Cleanup and replace ROOT types with basic types

* WARN: Disables misalignment when run under ROOT6.  Still
enabled under ROOT5.

* We don't even need a classdef here...

* We don't need to access most of this interactively... so comment out pragmas...

* Responsible for setting up the AgMLExtensions in the TGeoVolumes,
which connect user codes and volume IDs from the AgML source code
to the TGeo geometry model.

* Back to standard StMessMgr... ENDL-->endm

* Back to standard StMessMgr... ENDL-->endm

* The sanity check on the main branch is newer than the geant4star version.

* Expose pi to geoemtry modules

* Cleanup linkdef

* There is a pretty strong guarentee that table entries are not null
here b/c AgML would not be filling / creating null entries.  Plus...
it should be safe to delete a nullptr...

* Cautionary comments for future code maintainers...

* Rid ourselves of that commented-out classdef...

* Update StarVMC/StarAgmlLib/AgMath.h

* Apply suggestions from code review

* Fix memory leak in AgMLExtension

* Let's confirm the problem with ROOT6 first and then fix

https://github.com/star-bnl/star-sw/pull/308/files#r882223665

Co-authored-by: Dmitri Smirnov <[email protected]>

Author: klendathu2k

StarVMC/StarAgmlLib/AgBlock.cxx

@@ -5,12 +5,12 @@ ClassImp(AgBlock);
 #include <iostream>
 #include "AgMaterial.h"
 #include "AgShape.h"
-//#include "AgMedium.h"
 #include "AgModule.h"
 #include "AgPlacement.h"
+#include "AgMath.h" 
 
 #include "StarAgmlStacker.h"
-#include "StMessMgr.h"
+#include "StMessMgr.h"                                                                                                                  
 
 // Setup static members
 AgBlock                      *AgBlock::mCurrent=0;
@@ -51,7 +51,8 @@ AgBlock::AgBlock(const Char_t *name, const Char_t *title)
     mNicknames(),
     mGroups(),
     mMakeAssembly(0),
-    mNumberBranches(0)
+    mNumberBranches(0),
+    mGstpar()
 {
   addNickname(name); // ensure that the vector is populated
   _valid = true;
@@ -210,13 +211,3 @@ void AgBlock::List(Option_t *opts)
 
 }
 
-// ---------------------------------------------------------------------------------------------------
-// --
-//Bool_t AgBlock::isSame( const AgAttribute &attr )
-//{
-//  Bool_t same = true;
-//  same &= attr("serial") == _attribute("serial");
-//  return same;
-//
-//}
-

StarVMC/StarAgmlLib/AgBlock.h

@@ -14,7 +14,6 @@ class AgMixture;
 class AgVolume;
 class AgModule;
 #include "AgPlacement.h"
-#include "AgMath.h"
 #include "AgAttribute.h"
 #include "AgCreate.h"
 
@@ -154,6 +153,9 @@ class AgBlock : public TNamed
   int numberBranches() const { return mNumberBranches; }
   void branch(){ mNumberBranches++; }
 
+  void AddCut( TString cut, double value ){ mGstpar[cut] = value; }
+  std::map<TString, double>& GetCuts() { return mGstpar; }
+
  private:
  protected:
 
@@ -186,6 +188,8 @@ class AgBlock : public TNamed
   /// Counts the number of volume branchings
   int mNumberBranches;
 
+  std::map<TString, double> mGstpar; // GSTPAR cuts and parameters
+
   friend class _AgBlockDummy;
 
   ClassDef(AgBlock,0);

StarVMC/StarAgmlLib/AgMLExtension.cxx

@@ -1,6 +1,15 @@
 #include "AgMLExtension.h"
 
-AgMLExtension::AgMLExtension() : TGeoRCExtension(), mModuleName("none"), mFamilyName("none"), mVolumeName("none"), mSensitive(0), mTracking(0)
+AgMLExtension::AgMLExtension() : TGeoRCExtension(), 
+				 mModuleName("none"), 
+				 mFamilyName("none"), 
+				 mVolumeName("none"), 
+				 mSensitive(0), 
+				 mTracking(0),
+				 mBranchings(0),
+				 mVolumeId( new AgMLVolumeId ),
+				 mHitScoring(),
+				 mGstpar()
 {
   Grab();
 }

StarVMC/StarAgmlLib/AgMLExtension.h

@@ -3,23 +3,63 @@
 
 #include <TGeoExtension.h>
 #include <TString.h>
+#include <TMath.h>
+#include <vector>
+#include <map>
+
+// Volume ID functor (for sensitive volumes)
+class AgMLVolumeId {
+public:
+  virtual int id( int* numbv ) const { return 0; }
+};
+
+// User-defined hit scoring functor (for sensitive volumes)
+class AgMLScoring {
+public:
+  virtual float hit() const { return 0; }
+};
 
 class AgMLExtension : public TGeoRCExtension {
 
 public:
 
   AgMLExtension();
-  virtual ~AgMLExtension(){ /* nada */ };
+  virtual ~AgMLExtension() { delete mVolumeId; mVolumeId = nullptr; }
 
   void SetModuleName( const char* name ){ mModuleName = name; }
   void SetFamilyName( const char* name ){ mFamilyName = name; }
   void SetVolumeName( const char* name ){ mVolumeName = name; }
   void SetSensitive( bool flag ) { mSensitive = flag; }
-  void SetTracking( short track ) { mTracking = track; }
+  void SetTracking( short track ) { 
+    if ( mVolumeName == "CAVE" ) {
+      mTracking = 2;
+    }
+    else { 
+      mTracking = track; 
+    };
+  }
+
+  void SetBranchings( int b ) { mBranchings=b; }  
+
+  void SetVolumeIdentifier( AgMLVolumeId* identifier ){ delete mVolumeId; mVolumeId = identifier; }
+  void AddHitScoring( AgMLScoring* sc ){ mHitScoring.push_back( sc ); }
+
 
   TString GetModuleName(){ return mModuleName; }
   TString GetFamilyName(){ return mFamilyName; }
   TString GetVolumeName(){ return mVolumeName; }
+
+  int GetVolumeId( int* numbv ){ return mVolumeId ? mVolumeId->id( numbv ) : 0; }
+
+  bool GetSensitive() { return mSensitive; }
+  short GetTracking() { return mTracking; }
+  int GetBranchings(){ return mBranchings; }
+
+  std::vector<AgMLScoring*> GetUserHits(){ return mHitScoring; }
+
+  void AddCut( TString cut, double value ){ mGstpar[cut] = value; }
+  std::map<TString,double>& GetCuts() { return mGstpar; }
+  
 
 private:
 protected:
@@ -30,6 +70,12 @@ class AgMLExtension : public TGeoRCExtension {
 
   bool    mSensitive;  // volume sensitivity
   short   mTracking;   // 0=blackhole, 1=calorimeter, 2=tracking region
+  int     mBranchings; // number of branchings (placed family members)
+
+  AgMLVolumeId* mVolumeId; // Functor to calculate volume ID given reduced numbering scheme
+
+  std::vector<AgMLScoring*> mHitScoring; // Vector of functors for hit scoring
+  std::map<TString, double> mGstpar;     // GSTPAR tracking cuts for this volume
 
   ClassDef(AgMLExtension,0);
 

StarVMC/StarAgmlLib/AgMLStructure.h

@@ -220,7 +220,12 @@ class AgMLStructure : public AgMLStructureBase
 public:
 
   AgMLStructure( const char *_name ) : AgMLStructureBase(), name(_name) { };
-
+ ~AgMLStructure() {   
+   for ( auto t : table ) { // cleanup the table
+     delete (t);
+   }
+ };
+  
   void fill() {
     table.push_back( new T( current ) );
     int index = table.size() - 1; // current FILL index
@@ -261,30 +266,6 @@ class AgMLStructure : public AgMLStructureBase
     return false; 
   };
 
-
-  // //
-  // // Could simplify selextion by using c++11 lambdas...
-  // //
-  // //   HEXG.use( []( hexg_t h ){ return h.type == 1; } )
-  // //
-  // template<class Struct>
-  // bool use2( std::function<bool(Struct)> predicate ) {
-  //   auto x = std::find_if( table.begin(), table.end(), predicate );
-  //   if ( x == table.end() ) return false;
-
-  //   if ( *x ) { 
-  //     current = *(*x);
-  //     return true;
-  //   };
-
-  //   return false;
-  // };
-
-
-
-
-
-
   /// Cast to a refence to the current version
   operator T&(){ return current; }
 

StarVMC/StarAgmlLib/AgMaterial.cxx

@@ -5,23 +5,13 @@ ClassImp(AgMaterial);
 #include <assert.h>
 #include "AgBlock.h"
 #include "AgModule.h"
-#include "StMessMgr.h"
+#include "StMessMgr.h"                                                                                                                  
 
-//std::vector< TString >            AgMaterial::mParameterList;
 std::map< TString, AgMaterial * > AgMaterial::mMaterialTable;
 std::vector< TString >            AgMaterial::mMaterialList;
 
 struct _MaterialDummy { // Initializes material database
   _MaterialDummy(){
-    /* lifted
-    AgMaterial::mParameterList.push_back("a");
-    AgMaterial::mParameterList.push_back("z");
-    AgMaterial::mParameterList.push_back("dens");
-    AgMaterial::mParameterList.push_back("radl");
-    AgMaterial::mParameterList.push_back("absl");
-    AgMaterial::mParameterList.push_back("isvol");
-    AgMaterial::mParameterList.push_back("nelem");
-    */
 
     AgMaterial &h = AgMaterial::Get("Hydrogen");
     h.par("a")=1.010; h.par("z")=1.000; h.par("dens")=0.071; h.par("radl")=0.865E+03; h.par("absl")=0.790E+03; h.lock();
@@ -105,8 +95,6 @@ struct _MaterialDummy { // Initializes material database
     si.par("absl")=0.455E+02;    
     si.lock();
 
-//        18 LIQUID_ARGON            39.950    18.000     1.400 0.140E+02 0.837E+02   1
-
     AgMaterial &ar = AgMaterial::Get("Argon_gas");
     ar.par("a")=39.950;
     ar.par("z")=18.000;
@@ -150,8 +138,6 @@ struct _MaterialDummy { // Initializes material database
     mylar.Component("O",16.0,8.0,0.333);
     mylar.lock();
 
-    //$$$    AgMaterial::List();
-
   };
 
 } _material_dummy;
@@ -187,23 +173,6 @@ Double_t &AgMaterial::par( const Char_t *name ){
   return mParameters[name]; 
 }
 
-#if 0
-Bool_t AgMaterial::isSet( const Char_t *par ) const
-{
-  TString key=par;
-  return ( mParameters.find(key) != mParameters.end() );
-}
-
-Bool_t AgMaterial::hasPar(const Char_t *par ) const
-{
-  const char* key=par;
-  std::vector<std::string> parlist = mParameterList;
-  for ( UInt_t i=0;i<parlist.size();i++ )
-    if ( parlist[i] == key ) return true;
-  return false;
-}
-#endif
-
 Bool_t AgMaterial::hasComponent( const Char_t *comp )const
 {
 	if ( mType == AgMaterial::kMaterial )
@@ -216,7 +185,6 @@ void AgMaterial::Print( Option_t *opts ) const
 {
   TString name   = GetName();
   TString output = "";
-  //LOG_INFO << "["<< ((mLock)?"-":"+") << "] " << Form("%20s:",name.Data()) << " ";
   output = TString("[") + TString((mLock)?"-":"+") + TString("] ") + Form("%20s:",name.Data()) + TString(" ");
   std::map<std::string, Double_t > mypar=mParameters;
   for ( UInt_t j=0;j<mParameterList.size();j++ )
@@ -257,8 +225,6 @@ AgMaterial &AgMaterial::Get( const Char_t *name )
   TString modname="None";
   if ( module ) modname = module->GetName();
 
-  //  std::cout << "Info in <AgMaterial::Get(name)>: module="<<modname.Data()<<" material="<<name<<std::endl;
-
   //////////////////////////////////////////////////////////////
   //
   // If there is no module, probably being intialized above.
@@ -322,23 +288,16 @@ AgMaterial &AgMaterial::Get( const Char_t *name )
   //
   ////////////////////////////////////////////////////////////// 
 
-  //std::cout << "Returning new material " << matname.Data() << std::endl;
-
   material = new AgMaterial(matname);
 
   if ( AgBlock::previous() )                    // Inherit from creating block
     material->Inherit( AgBlock::previous() ); 
   else if ( AgBlock::module() )                 // Or from the current material loaded in module    
     material->Inherit( AgBlock::module() );
 
-  //std::cout << "++ After inheritance from mother block" << std::endl;
-
   if ( AgBlock::active() )                      // And finally the active block
     material->Inherit( AgBlock::active() );   
 
-  //std::cout << "++ After inheritance from active block" << std::endl;
-  //material->Print();
-
   material->SetName(matname);
   mMaterialTable[matname]=material;
   mMaterialList.push_back(matname); 
@@ -360,9 +319,6 @@ AgMaterial AgMaterial::CopyMaterial( const Char_t *name )
   if ( module  ) modname = module->GetName();
   if ( current ) volname = current->GetName();
 
-
-  //  std::cout << "Info in <AgMaterial::Get(name)>: module="<<modname.Data()<<" material="<<name<<std::endl;
-
   //////////////////////////////////////////////////////////////
   //
   // If there is no module, probably being intialized above.

StarVMC/StarAgmlLib/AgMaterial.h

@@ -7,7 +7,6 @@
 #include <map>
 #include <vector>
 
-//#include "AgParameterList.h"
 #include <StarVMC/StarAgmlUtil/AgParameterList.h>
 class AgMaterial;
 class AgBlock;
@@ -23,17 +22,9 @@ class AgMaterial : public TNamed, public AgParameterList<double>
 
   enum { kUnknown=0, kMaterial, kMixture, kCompound };
 
-
   /// Returns a reference to the named parameter.
   Double_t &par( const Char_t *name );
 
-#if 0 // lifted
-  /// Tests whether the named parameter is set for this shape
-  Bool_t isSet( const Char_t *par ) const;  
-  /// Tests whether the named parameter is valid for this shape
-  Bool_t hasPar( const Char_t *par ) const;
-#endif
-
   /// Tests wether the mixture has the named component.  Returns false
   /// if this is a material.
   Bool_t hasComponent( const Char_t *comp ) const;
@@ -75,7 +66,6 @@ class AgMaterial : public TNamed, public AgParameterList<double>
 
   /// Returns (by reference) the name, A, Z and Weight of each component
   void Component( Int_t ic, TString &name, Double_t &a, Double_t &z, Double_t &weight );
-
 
   Bool_t locked(){ return mLock; }
 
@@ -85,6 +75,13 @@ class AgMaterial : public TNamed, public AgParameterList<double>
   void   lock(){ mLock=true; }
 
   Double_t sumWeights();
+
+  struct MyComponent {
+	  Double_t a;
+	  Double_t z;
+	  Double_t w;
+  };
+
 
  private:
  protected:
@@ -98,12 +95,6 @@ class AgMaterial : public TNamed, public AgParameterList<double>
   std::vector< TString >  mC;
   std::vector< Double_t > mA, mZ, mW;
 
-  struct MyComponent {
-	  Double_t a;
-	  Double_t z;
-	  Double_t w;
-  };
-
   std::map< TString, MyComponent > mComponentTable; // Table storing already defined components
 
   static std::map< TString, AgMaterial * > mMaterialTable; // List of existing materials