MCApplication.cxx

Go to the documentation of this file.

//
// MCApplication
//
// MCApplication is a pure virtual base class to be used to define the required
// interface of derived vmc Application classes, e.g. PTSimApplication
// and GeoSwimmerApplication.
//
// Author:  S. Kasahara 10/06
//
// Notes: This class contains a subset of the methods defined in
//        $ROOTSYS/vmc/inc/TVirtualMCApplication class.
//

#include <TVirtualMC.h>
#include <TInterpreter.h>
#include <TRandom3.h>
#include <TROOT.h>
#include <TSystem.h>
#include <TDatabasePDG.h>

#include "Util/UtilPDG.h"
#include "Util/UtilMCFlag.h"
using namespace UtilMCFlag;
#include "Util/LoadMinosPDG.h"
#include "MCApplication/MCApp.h" // for process/cut flags
#include "MCApplication/MCApplication.h"
#include "MessageService/MsgService.h"
#include "BField/BField.h"
#include "GeoGeometry/GeoMedium.h"

ClassImp(MCApplication)

CVSID("$Id: MCApplication.cxx,v 1.21 2009/06/22 03:51:00 kasahara Exp $");

//_____________________________________________________________________________
const MCApplication& MCApplication::Instance() {
  // Access to singleton instance of MCApplication

  // fgInstance belongs to base class
  if ( !TVirtualMCApplication::Instance() ) {
    // Construct an Instance of this
    static MCApplication::Cleaner c;
    c.ClassIsUsed();
    new MCApplication("MCApplication","MINOS MC Application");
  }
  
  return dynamic_cast<const MCApplication&>
                     (*(TVirtualMCApplication::Instance()));
  
}
  
//_____________________________________________________________________________
MCApplication::MCApplication(const char *name, const char *title) 
  : TVirtualMCApplication(name,title), fMCAppUser(0),fMC(0),
    fUgliGeomHandle(0),fVldContext(),fBField(0),fBFldInZTestSave(0),fRandom(0),
    fMCConfigScript(""),fDecayConfigScript(""),fLogLevel(kInfo) {
  // Normal constructor.  Private.  Use MCApplication::Instance() to
  // construct Instance of MCApplication singleton.

  MSG("MCApp",Msg::kDebug) << "MCApplication normal ctor @ " << this << endl;

  // Initialize random number generator
  fRandom = new TRandom3();

  LoadMinosPDG(); // initialize MINOS specific particles
  
}


//_____________________________________________________________________________
MCApplication::~MCApplication() {
  // Destructor.  Delete all owned sub-objects

  MSG("MCApp",Msg::kDebug) << "MCApplication dtor @ " << this << endl;

  Clear(); // clear dynamically allocated memory

}

//_____________________________________________________________________________
void MCApplication::Clear(const Option_t* /* option */) {
  // Protected method to Clear dynamically allocated memory.  

  // fMCAppUser not owned.
  if ( fUgliGeomHandle ) delete fUgliGeomHandle; fUgliGeomHandle = 0;
  if ( fBField ) delete fBField; fBField = 0;
  if ( fRandom ) delete fRandom; fRandom = 0;
  if ( fMC ) delete fMC; fMC = 0;

}

//_____________________________________________________________________________
void MCApplication::SetUserApplication(MCAppUser* appuser) {
  // Set new user application to be called during particle tranport

  fMCAppUser = appuser;
  if ( fMC ) fMC -> SetStack(fMCAppUser->GetStack());
 
}
  
//_____________________________________________________________________________
void MCApplication::InitMC(const char* mcconfigscript, const VldContext& vldc,
                           const char* decayconfigscript) {
  // Initialize MC through execution of configuration script. Can only
  // be done once.   Also initialize geometry appropriate for input validity
  // context.  

  // Store vldcontext to be used in construction of geometry
  fVldContext = vldc;
  fMCConfigScript = mcconfigscript;
  fDecayConfigScript = decayconfigscript;
  fLogLevel = MsgService::Instance()->GetStream("MCApp")->GetLogLevel();
  
  // Initialize MC from configuration script
  if ( fMC ) { 
    MSG("MCApp",Msg::kWarning) << "MCApplication::InitMC called with script "
       << fMCConfigScript.c_str() << "\nbut MC already initialized. Ignored." 
       << endl;
    return;
  }
  if ( fMCConfigScript.empty() || fMCConfigScript == std::string("<null>") ) {
    MSG("MCApp",Msg::kFatal) << "MCApplication::InitMC called with null "
                             << " configuration script." << endl;
    abort();
  }
  
  // Print, load, and execute configuration script
  std::string mcConfigScriptPath = gSystem->Which(gROOT->GetMacroPath(),
                                   fMCConfigScript.c_str(),kReadPermission);
  if ( mcConfigScriptPath.empty() ) {
    MSG("MCApp",Msg::kError) << "MCApplication::InitMC configuration script "
                           << fMCConfigScript.c_str() << " not found in path "
                           << gROOT->GetMacroPath() << ". Abort." << endl;
    abort();
  }
    
  MSG("MCApp",Msg::kInfo) 
          << "MCApplication::InitMC initializing MC with configuration script "
          << mcConfigScriptPath.c_str() << ":" << endl;

  if ( fLogLevel <= Msg::kInfo ) {
    std::string printstr = ".!cat -n " + mcConfigScriptPath;
    gInterpreter->ProcessLineSynch(printstr.c_str());
  }
  gROOT->LoadMacro(mcConfigScriptPath.c_str());
  gInterpreter->ProcessLineSynch("Config()");
  
  if ( !gMC ) {
    MSG("MCApp",Msg::kFatal) << "Construction of MC implementation failed."
                             << endl;
    abort();
  }
  fMC = gMC;

  // Now begin initialization of MC

  // random number generator can be seeded in BeginEvent by user app
  fMC -> SetRandom(fRandom); 
  if ( fMCAppUser ) fMC -> SetStack(fMCAppUser->GetStack());
  // Init() will call MCApplication::AddParticles() to define additional
  // particles not defined in MC by default and 
  // MCApplication::ConstructGeometry to build geometry
  fMC -> Init(); 

  // Initialize parameters of swim media for continuous energy loss
  // This has to be done after Init() method has been called, but before
  // BuildPhysics();
  InitMedia();
  
  fMC -> BuildPhysics();
  
}

//_____________________________________________________________________________
bool MCApplication::InitSnarl(const VldContext& vldc) {
  // Initialize Snarl according to current VldContext.  

  bool isOkay = true;

  if ( !fBField || !fUgliGeomHandle ) {
    MSG("MCApp",Msg::kError) << "InitSnarl w/" << vldc
                             << " called before InitMC. Abort." 
                             << endl;
    abort();
  }
  
  // Reset BField to match snarl vldcontext
  fBField -> ResetVldContext(vldc);

  // Check that Geometry hasn't changed.  
  const VldRange& vldrange = fUgliGeomHandle -> GetVldRange();
  
  if ( !(vldrange.IsCompatible(vldc)) ) {
    MSG("MCApp",Msg::kError) << "InitSnarl w/vld " << vldc
                             << " is incompatible with VMC Geometry vld "
                             << vldrange << ". Abort." << endl;
    abort();
  }
  
  return isOkay;
  
}

//_____________________________________________________________________________
void MCApplication::FatalNoApp() const {
  // Issue fatal warning message and abort

  MSG("MCApp",Msg::kFatal) << "MCApplication has null fMCAppUser. abort."
                           << endl;
  abort();
  
}

//_____________________________________________________________________________
void MCApplication::ConstructGeometry() {
  // Construct geometry.  Called once during initialization by 
  // TVirtualMC::Init()

  MSG("MCApp",Msg::kDebug) << "MCApplication::ConstructGeometry for vldc "
                           << fVldContext.AsString() << endl;
  
  if ( fUgliGeomHandle ) { 
    MSG("MCApp",Msg::kFatal) 
    << "ConstructGeometry called but fUgliGeomHandle non-null. abort." << endl;
    abort();
  }
  if ( !fVldContext.IsValid() ) {
    MSG("MCApp",Msg::kFatal)
      << "ConstructGeometry but invalid fVldContext. abort." << endl;
    abort();
  }

  fUgliGeomHandle = new UgliGeomHandle(fVldContext);
  fBField = new BField(fVldContext);

  // Update gGeoManager to guarantee that it points to the right place before
  // returning
  UpdateGlobalGeoManager();

  // Notify VMC about ROOT geometry
  gMC -> SetRootGeometry();
 
  return;
  
}

//_____________________________________________________________________________
void MCApplication::AddParticles() {
  // Define particles beyond those defined by MC.  Called once during 
  // initialization by concrete MC.
  // Also define additional or redefine existing Geant3 decay modes.

  MSG("MCApp",Msg::kDebug) << "MCApplication::AddParticles" << endl;

  // Define ion types
  AddIons();

  // Define decay modes 
  AddDecayModes();

  return;
  
}

//_____________________________________________________________________________
void MCApplication::AddIons() {
  // Private method to define ions to be used in monte carlo simulation.
  // Called by AddParticles() during initialization.
  // Ions are defined according to ion particle definitions in
  // TDatabasePDG, as filled by LoadMinosPDG, using pdg-2006 particle
  // codes.

  MSG("MCApp",Msg::kDebug) << "MCApplication::AddIons" << endl;

  // Loop over particles defined in TDatabasePDG and define
  // those up to z = zmax to Monte Carlo
  const Int_t zmax = 28;

  const TDatabasePDG& dbpdg = *(TDatabasePDG::Instance());
  TIter next(dbpdg.ParticleList());
  TParticlePDG* particle = 0;
  Int_t ia,iz,ij;
  Double_t lifetime = 1.e20; // lifetime in sec, stable for all ions

  while ((particle = (TParticlePDG*)next())) {
    Int_t pdgcode = particle -> PdgCode();
    UtilPDG::ionscheme_t ionscheme = UtilPDG::getIonAZJ(pdgcode,ia,iz,ij);
    if ( ionscheme == UtilPDG::kPDG2006 && iz <= zmax ) {
      std::string pname = particle -> GetName();
      Double_t mass = particle -> Mass();
      Double_t charge = (particle -> Charge())/3.; //charge in units of |e|
      if ( gMC->IdFromPDG(pdgcode) <= 0 ) {
#if ROOT_VERSION_CODE >= ROOT_VERSION(5,22,0)
        gMC->DefineParticle(pdgcode,pname.c_str(),kPTIon,mass,charge,lifetime,
                           "nucleus", 0.0, 1, 1, 0, 1, 1, 0, 0, 1, kTRUE);
#else
        gMC->DefineParticle(pdgcode,pname.c_str(),kPTIon,mass,charge,lifetime);
#endif
        MSG("MCApp",Msg::kDebug) << "AddIon " << pdgcode << "/" 
           << pname.c_str() << " mass " << mass << " charge " << charge 
           << " lifetime " << lifetime << endl;
      }    
    }
  }

  // Add additional ion types with z > zmax to cover all those
  // defined by GMINOS.
  const Int_t nion = 15+12+1;
  // The first set of 15 ions are the subset of ions defined by 
  // Geant3 GPIONS .AND. TDatabasePDG.  (Ions defined in GPIONS but
  // not TDatabasePDG are not added and are commented out below.)   
  // The TDatabasePDG definitions of this set of ions is used in place
  // of that defined in GPIONS for consistency with the rest of the
  // ion definitions.
  // The second set of 12 ions are those defined by GMINOS lbl_part.
  // The third set of ions are those defined because they've popped
  // up during simulation.
  Int_t izia[nion][2] = {{29,63},  // Cu63
                         {30,64},  // Zn64
                         {32,74},  // Ge74
                         {34,80},  // Se80
                         {36,84},  // Kr84 
                         //  {38,88},  // Sr88
                         //  {40,90},  // Zr90
                         //  {42,98},  // Mo98
                         {46,106}, // Pd106
                         {48,114}, // Cd114
                         {50,120}, // Sn120 
                         {54,132}, // Xe132
                         {56,138}, // Ba138
                         //  {58,140}, // Ce140
                         //  {62,152}, // Sm152
                         //  {66,164}, // Dy164
                         //  {70,174}, // Yb174
                         {74,184}, // W184
                         {78,194}, // Pt194
                         {79,197}, // Au197
                         {80,202}, // Hg202
                         {82,208}, // Pb208
                         //  {92,238}, // U238
                         {31,69},  // Ga69
                         {33,75},  // As75
                         {35,79},  // Br79
                         {47,107}, // Ag107
                         {49,115}, // In115
                         {51,121}, // Sb121
                         {52,130}, // Te130
                         {53,127}, // I127
                         {55,133}, // Cs133
                         {73,181}, // Ta181
                         {77,193}, // Ir193
                         {81,205}, // Tl205
                         {50,115}, // Sn115
                                 };
  
  for ( int iion = 0; iion < nion; iion++ ) {
    Int_t pdgcode = UtilPDG::makeIonPDG(izia[iion][1],izia[iion][0],ij=0,
                                         UtilPDG::kPDG2006);
    TParticlePDG* particle = dbpdg.GetParticle(pdgcode);
    if ( particle ) {
      std::string pname = particle -> GetName();
      Double_t mass = particle -> Mass();
      Double_t charge = (particle -> Charge())/3.; //charge in units of |e|
      if ( gMC->IdFromPDG(pdgcode) <= 0 ) {
#if ROOT_VERSION_CODE >= ROOT_VERSION(5,22,0)
        gMC->DefineParticle(pdgcode,pname.c_str(),kPTIon,mass,charge,lifetime,
                           "nucleus", 0.0, 1, 1, 0, 1, 1, 0, 0, 1, kTRUE);
#else
        gMC->DefineParticle(pdgcode,pname.c_str(),kPTIon,mass,charge,lifetime);
#endif
        MSG("MCApp",Msg::kDebug) << "AddIon " << pdgcode << "/" 
           << pname.c_str() << " mass " << mass << " charge " << charge 
           << " lifetime " << lifetime << endl;
      }    
    }
    else {
      MSG("MCApp",Msg::kWarning) << "MCApplication::AddIons. Ion " 
                                 << pdgcode << " not present "
                                 << "in TDatabasePDG. Will not be defined!"
                                 << endl;
    }
  }
  
  return;
  
}

//_____________________________________________________________________________
void MCApplication::AddDecayModes() {
  // Add decay modes beyond or overriding those defined by the 
  // concrete VMC.
  // This is done through the activation of an external decayer, e.g.
  // TPythia6Decayer, through an external script. An example of such
  // a script is PTSim_DecayConfig.C.
  
  if ( fDecayConfigScript.empty() 
       || fDecayConfigScript == std::string("<null>") ) {
    MSG("MCApp",Msg::kDebug) << "MCApplication::AddDecayModes. No "
                         << "external decay script defined." << endl;
    return;
  }
  
  // Print, load, and execute configuration script
  std::string decayConfigScriptPath = gSystem->Which(gROOT->GetMacroPath(),
                               fDecayConfigScript.c_str(),kReadPermission);
  if ( decayConfigScriptPath.empty() ) {
    MSG("MCApp",Msg::kError) << "MCApplication::InitMC configuration script "
                        << fDecayConfigScript.c_str() << " not found in path "
                        << gROOT->GetMacroPath() << ". Abort." << endl;
    abort();
  }
    
  // Configure decayer using an external configuration script
  MSG("MCApp",Msg::kInfo) 
   << "MCApplication::AddDecayModes initializing decayer with script "
   << decayConfigScriptPath.c_str() << ":" << endl;

  if ( fLogLevel <= Msg::kInfo ) {
    std::string printstr = ".!cat -n " + decayConfigScriptPath;
    gInterpreter->ProcessLineSynch(printstr.c_str());
  }
  gROOT->LoadMacro(decayConfigScriptPath.c_str());
  gInterpreter->ProcessLineSynch("DecayConfig()");

  return;
  
}

//_____________________________________________________________________________
void MCApplication::InitGeometry() {
  // Initialize geometry. Called by TVirtualMC::Init().  

  return;
  
}

//_____________________________________________________________________________
void MCApplication::InitMedia() {
  // Private method to initialize parameters of media according to
  // user configured cut & process values if specified.

  UpdateGlobalGeoManager(); 
  
  TIter nextmed(gGeoManager->GetListOfMedia());
  TGeoMedium* med = 0;
  while ( ( med = (TGeoMedium*)nextmed() ) ) {
    GeoMedium* geomed = dynamic_cast<GeoMedium*>(med);
    if ( !geomed ) continue;
    std::string medName = geomed -> GetName();
    const GeoMedium::CutMap& cutmap = geomed -> GetCutMap();
    for ( GeoMedium::CutMapConstItr citr  = cutmap.begin(); 
                                    citr != cutmap.end(); citr++ ) {
      MSG("Geo",Msg::kDebug) << "Medium " << medName.c_str()
                             << " " << UtilMCFlag::AsString(citr->first)
                             << " set to " << citr->second 
                             << " overriding default." << endl;
      fMC -> Gstpar(geomed->GetId(),UtilMCFlag::AsString(citr->first),
                                    citr->second);
    }
    const GeoMedium::ProcessMap& processmap = geomed -> GetProcessMap();
    for ( GeoMedium::ProcessMapConstItr citr = processmap.begin(); 
                                        citr!= processmap.end(); citr++ ) {
      MSG("Geo",Msg::kDebug) << "Medium " << medName.c_str()
                             << " " << UtilMCFlag::AsString(citr->first)
                             << " set to " << citr->second 
                             << " overriding default." << endl;
      fMC -> Gstpar(geomed->GetId(),UtilMCFlag::AsString(citr->first),
                                    citr->second);
    }
  }

}

//_____________________________________________________________________________
void MCApplication::Field(const Double_t* x, Double_t* b) const {
  // Determine magnetic field b[3] at global position x[3]
  //
  // Units of BField map are meters for input position and Tesla for field
  // Units of MC application are cm for position and kiloGauss for field.
  //

  // Assume null
  b[0] = 0; b[1] = 0; b[2] = 0;
  
  if ( !fBField ) {
    MSG("MCApp",Msg::kWarning) << "MCApplication::Field but Null fBField!"
                               << endl;
    return;
  }

  // Insurance that current medium has field activated, although in principle 
  // this shouldn't be necessary 
  TGeoNode* node = gGeoManager->GetCurrentNode();
  if ( node -> GetVolume() -> GetMedium() -> GetParam(1) <= 0 ) return;
  
  TVector3 bxyz(0,0,0);
  TVector3 xyz(x[0]*Munits::cm,x[1]*Munits::cm,x[2]*Munits::cm);

  bxyz = fBField -> GetBField(xyz);

  // Convert Tesla to kGauss
  for ( int i = 0; i < 3; i++ ) b[i] = bxyz[i]*10.;

  UpdateGlobalGeoManager(); // shouldn't be necessary, but caution rules

  return;
 
  
}

//_____________________________________________________________________________
void MCApplication::BeginEvent() {
  // Action to be taken at beginning of an event
  //

  UpdateGlobalGeoManager(); // in case gGeoManager has been altered

  // Deactivate BField test to determine if test pt is within steel plane
  // when using VMC.  Store current state of test, to be restored
  // post particle transport in FinishEvent()
  fBFldInZTestSave = fBField -> GetRequireInZTest();
  //fBField -> SetRequireInZTest(0);
  
  if ( !fMCAppUser ) FatalNoApp();
  fMCAppUser -> BeginEvent();

  return;
  
}

//_____________________________________________________________________________
void MCApplication::FinishEvent() {
  // Action to be taken at end of an event

  UpdateGlobalGeoManager(); // in case gGeoManager has been altered

  // Restore BField test to status before particle transport
  fBField -> SetRequireInZTest(fBFldInZTestSave);

  if ( !fMCAppUser ) FatalNoApp();
  fMCAppUser -> FinishEvent();
  
  return;
  
}

  
  


  

---