PTSimModule.cxx

Go to the documentation of this file.

//
// PTSimModule
//
// PTSimModule is a JobControl Module for use with simulating particle
// transport through the MINOS detectors.
//
// S. Kasahara 05/04

#include <cassert>
//#include "/usr/include/google/profiler.h"

#include <TClonesArray.h>
#include <TVirtualMC.h>
#include <TMCProcess.h>
#include <TPDGCode.h>
#include <TDatabasePDG.h>

#include "MessageService/MsgService.h"
#include "JobControl/JobCommand.h"
#include "JobControl/JobCModuleRegistry.h"
#include "ParticleTransportSim/PTSim.h"
#include "ParticleTransportSim/PTSimModule.h"
#include "ParticleTransportSim/PTSimApplication.h"
#include "ParticleTransportSim/PTSimStack.h"
#include "MCApplication/MCApplication.h"
#include "Util/LoadMinosPDG.h"
#include "Util/UtilString.h"
#include "Record/SimSnarlRecord.h"
#include "MinosObjectMap/MomNavigator.h"

ClassImp(PTSimModule)

//   Definition of static data members
//   *********************************

CVSID("$Id: PTSimModule.cxx,v 1.56 2009/06/22 03:51:01 kasahara Exp $");

JOBMODULE(PTSimModule, "PTSimModule",
          "A module for running the particle transport simulation.");

//_____________________________________________________________________________
PTSimModule::PTSimModule(): fMCConfig("<null>"),fDecayConfig("<null>"),
                            fMCAppUser(0),fRandomSeed(0) {
  // Default constructor

  LoadMinosPDG(); // Load PDG database w/MINOS specific particles
  
  fStopwatch.Reset();
  fStopwatch.Stop();
  

  // The application must be instantiated before the mc implementation
  fMCAppUser = new PTSimApplication("PTSimApplication","Minos MC Application");
  if ( !fMCAppUser ) {
    MSG("PTSim",Msg::kFatal) << "No MC application instance." << endl;
    abort();
  }


}

//_____________________________________________________________________________
PTSimModule::~PTSimModule() {
   // Destructor

  MSG("PTSim",Msg::kDebug) << "PTSimModule dtor @ " << this << endl;
  
  if ( fMCAppUser ) delete fMCAppUser; fMCAppUser = 0;

  // Clear contents of MCApplication because MCApplication singleton
  // destruction is out of order with G4 singleton destruction
  MCApplication& mcapp = const_cast<MCApplication&>(MCApplication::Instance());
  mcapp.Clear();
  
}

//_____________________________________________________________________________
const Registry& PTSimModule::DefaultConfig() const {
  //
  // Purpose: Method to return default configuration.
  // 
  // Arguments: none.
  //
  // Return: Registry containing default configuration
  //

  MSG("PTSim",Msg::kDebug) << "Loading default config" << endl;

  static Registry r; 
  std::string name = this->JobCModule::GetName();
  name += ".config.default";
  r.SetName(name.c_str());

  r.UnLockValues();
  r.Set("MCConfig","<null>");
  r.Set("DecayConfig","<null>");  // or PTSim_DecayConfig.C
  r.Set("StdHepSelectMask",PTSim::kMomentum);
  // stdhep select mask default is to save
  // all tracks above momentum threshold (StdHepThr) 
  r.Set("StdHepThr",0.15);// default p(GeV/c) thresh applied to stdhep trks
  // if user has selected mode to save secondaries resulting in energy
  // deposition hit, StdHepHitThr determines threshold above which to
  // store
  r.Set("StdHepHitThr",0.001);// default E(GeV) hit threshold
  // flag used to determine if selected secondary siblings are to be saved 
  // (0 => false, 1 => true, 2 (default) => true and save siblings of direct 
  // ancestors as well if StdHepSaveAncestor true)
  r.Set("StdHepSaveSibling",2);// default is to save selected secondary sibs
                               // and sibs of ancestors
  r.Set("StdHepSaveAncestor",1);// default is to save selected secondary ances
  r.Set("RandomSeed",fRandomSeed);  // optional random seed supplied by user
  r.Set("RkVLevel",0);  // default is to impose no rock veto cuts
  r.Set("RkVMinDist",2.); // default minimum distance at which veto cut applies
  // safety factor multiplied in front of E/(dE/dX)_min when determining max 
  // distance particle will travel 
  r.Set("RkVFactor",1.25); 
  r.Set("UseLinearStep",false); // bad idea to set true, but useful for validation against gminos

  r.LockValues();

  // Set defaults for storing secondaries
  // Always store secondary muons, regardless of production mechanism or 
  // momentum
  PTSimStack* stack = fMCAppUser -> GetStack();
  if ( !stack ) {
    MSG("PTSim",Msg::kFatal) << "No stack!" << endl;
    abort();
  }

  // Always store muons at production
  stack -> SetStdHepThrByType(0.,kPNoProcess,-13);
  stack -> SetStdHepThrByType(0.,kPNoProcess,+13);
  // Always store secondaries produced by muon decay, regardless of type
  stack -> SetStdHepThrByType(0.,kPDecay,kRootino,-13);
  stack -> SetStdHepThrByType(0.,kPDecay,kRootino,+13);


  return r;
}


//_____________________________________________________________________________
void PTSimModule::Config(const Registry& r) {
  //
  // Purpose: Configure the module given a registry.
  //
  // Arguments: Registry to use to configure the module.
  //
  // Return: none.
  //

  MSG("PTSim",Msg::kDebug) << "Config PTSimModule with r=" << r << endl;
  
  Int_t tmpi;
  Double_t tmpd;
  const Char_t* tmpc;
  if ( r.Get("MCConfig",tmpc) ) fMCConfig = tmpc;
  if ( r.Get("DecayConfig",tmpc) ) fDecayConfig = tmpc;

  PTSimStack* stack = fMCAppUser -> GetStack();
  if ( r.Get("StdHepThr",tmpd) ) stack -> SetStdHepThr(tmpd);
  if ( r.Get("StdHepHitThr",tmpd) ) stack -> SetStdHepHitThr(tmpd);
  if ( r.Get("StdHepSelectMask",tmpi) ) stack -> SetStdHepSelectMask(tmpi);
  if ( r.Get("StdHepSaveSibling",tmpi) ) stack -> SetStdHepSaveSibling(tmpi);
  if ( r.Get("StdHepSaveAncestor",tmpi) ) stack -> SetStdHepSaveAncestor(tmpi);
  if ( r.Get("RandomSeed",tmpi) ) fRandomSeed = tmpi;
  if ( r.Get("RkVLevel",tmpi) ) fMCAppUser -> SetRkVLevel(tmpi);
  if ( r.Get("RkVMinDist",tmpd) ) fMCAppUser -> SetRkVMinDistInCM(tmpd*100.);
  if ( r.Get("RkVFactor",tmpd) ) fMCAppUser -> SetRkVFactor(tmpd);
  if ( r.Get("UseLinearStep",tmpi) ) fMCAppUser -> UseLinearStep(tmpi);
  
}

//_____________________________________________________________________________
void PTSimModule::BeginJob() {
  //
  //  Purpose:  Called once at the beginning of the job to
  //            handle job module initialization.
  //
  //  Arguments: none.
  //  
  //  Return: status.
  // 

   MSG("PTSim",Msg::kDebug) << "PTSimModule::BeginJob" << endl;

   // mc implementation initialization is saved until first event
   // since vldc is required for geometry construction
   
}

//_____________________________________________________________________________
void PTSimModule::EndJob() {
  //
  //  Purpose:  Called once at the end of the job to handle job
  //            module validation.
  //
  //  Arguments: none.
  //  
  //  Return: status.
  // 

   MSG("PTSim",Msg::kDebug) << "PTSimModule::EndJob" << endl;

   //ProfilerStop();
   
   fStopwatch.Stop();
   Int_t nSnarl = 0;
   Int_t nEvent = 0;
   if ( fMCAppUser ) {
     nSnarl = fMCAppUser -> GetNSnarls();
     nEvent = fMCAppUser -> GetNEvents();
   }
   
   MSG("PTSim",Msg::kInfo) << "PTSimModule::EndJob, Time(sec), Real "
                           << fStopwatch.RealTime() << ", CPU "
                           << fStopwatch.CpuTime() << " # snarls/events " 
                           << nSnarl << "/" << nEvent <<  endl;
   
}

//_____________________________________________________________________________
void PTSimModule::HandleCommand(JobCommand* cmd) {
  // Used to handle commands that don't work Config
  //
  
  if ( !cmd->HaveCmd() ) return;
  string s = cmd->PopCmd();
  if ( s == "StdHepThrByType" ) {
    UtilIstHEP::EProdMethod prodmethod = UtilIstHEP::kMNoProcess;
    Int_t pdgId = kRootino;
    Int_t parentId = kRootino;
    double thresh = 0.;
    if ( cmd->HaveOpt() ) {
      thresh = cmd->PopFloatOpt();

      if ( cmd->HaveOpt() ) {
        int prodmethodint = cmd -> PopIntOpt();
        if ( prodmethodint < 0 || prodmethodint >= UtilIstHEP::kNProdMethod ){
          MSG("PTSim",Msg::kWarning) 
          << "StdHepThrByType called w/invalid UtilIstHEP::EProdMethod "
          << prodmethodint << ". Ignored." << endl;
          return;
        }
        else prodmethod = (UtilIstHEP::EProdMethod)(prodmethodint);

        if ( cmd->HaveOpt() ) {
          pdgId = cmd -> PopIntOpt();
          if ( cmd->HaveOpt() ) parentId = cmd -> PopIntOpt();
        }
      }
      else {
        MSG("PTSim",Msg::kWarning)
          << "StdHepThrByType called w/no production method. Ignored." 
          << endl;
        return;
      }  
    }
    else {
      MSG("PTSim",Msg::kWarning)
        << "StdHepThrByType called w/no momentum threshold. Ignored."
        << endl;
      return;
    }
    
    // Particle Name
    std::string particlename = " for particles of type ";
    if ( pdgId != kRootino ) {
      particlename += UtilString::ToString<int>(pdgId) + "/";
      const TDatabasePDG& dbpdg = *(TDatabasePDG::Instance());
      particlename += "???";
      if ( dbpdg.GetParticle(pdgId) ) particlename 
                       += dbpdg.GetParticle(pdgId)->GetName();
    }
    else particlename = "";

    std::string parentname = " for particles w/parent type ";
    if ( !particlename.empty() ) parentname = " and parent type ";
    if ( parentId != kRootino ) {
      parentname += UtilString::ToString<int>(parentId) + "/";
      const TDatabasePDG& dbpdg = *(TDatabasePDG::Instance());
      parentname += "???";
      if ( dbpdg.GetParticle(parentId) ) parentname 
                       += dbpdg.GetParticle(parentId)->GetName();
    }
    else parentname = "";
    
    if ( prodmethod == UtilIstHEP::kMNoProcess ) {
      MSG("PTSim",Msg::kInfo) << "PTSimModule StdHepThrByType " 
                              << "configured to " << thresh << "(GeV/c)"
                              << particlename.c_str()
                              << parentname.c_str()
                              << endl;
    }
    else {
      std::string prodmethodname = " for particles produced by ";
      if ( !parentname.empty() || !particlename.empty() ) 
        prodmethodname = " and produced by ";
      MSG("PTSim",Msg::kInfo) << "PTSimModule StdHepThrByType " 
                              << "configured to " << thresh << "(GeV/c)"
                              << particlename.c_str()
                              << parentname.c_str()
                              << prodmethodname.c_str() << prodmethod 
                              << "/" << UtilIstHEP::AsString(prodmethod) 
                              << endl;
    }
    
    PTSimStack* stack = fMCAppUser -> GetStack();
    TMCProcess tmcprocess = UtilIstHEP::GetTMCProcess(prodmethod);
    stack -> SetStdHepThrByType(thresh,tmcprocess,pdgId,parentId);
    return;
  }
  else if ( s == "StdHepSave" ) {

    std::vector<int> snarllist;
    MsgStream& msgDebug = MSGSTREAM("PTSim",Msg::kDebug);
    msgDebug << "StdHepSave ";

    bool isrange = false;
    Int_t nsnarl = 0;
    while ( cmd->HaveOpt() ) {
      Int_t snarlno = cmd->PopIntOpt();
      msgDebug << snarlno << " ";
      snarllist.push_back(TMath::Abs(snarlno));
      nsnarl++;
      if ( snarlno < 0 ) isrange = true;
    }
    msgDebug << endl;
    if ( isrange && nsnarl != 2 ) {
      MSG("PTSim",Msg::kWarning)
        << "StdHepSave by range requires two arguments. Ignored." << endl;
      snarllist.clear();
      isrange = false;
    }
    
    PTSimStack* stack = fMCAppUser -> GetStack();
    stack -> SetStdHepSave(snarllist,isrange);

    return;
  }

  // Errors fall through to here
  MSG("PTSim",Msg::kError) << "Illegal command: " << s.c_str() << endl;
  
}

//_____________________________________________________________________________
JobCResult PTSimModule::Reco(MomNavigator *mom) {
  //
  //  Purpose:  Run the requested mc implementation.
  //
  //  Arguments: mom.
  //  
  //  Return: status.
  // 

  JobCResult result(JobCResult::kPassed);  
  MSG("PTSim",Msg::kDebug) << "PTSimModule::Reco" << endl;

  // Check that mom exists.
  assert(mom);

  SimSnarlRecord* record 
  = dynamic_cast<SimSnarlRecord*>(mom->GetFragment("SimSnarlRecord"));
  if ( !record ) {
    MSG("PTSim",Msg::kError) << "No SimSnarlRecord in Mom." << endl;
    result.SetWarning().SetFailed();
    return result;
  }
  
  RecJobHistory& jobhist 
              = const_cast<RecJobHistory&>(record->GetJobHistory());
  jobhist.CreateJobRecord(RecJobHistory::kPTSim);
  
  const SimSnarlHeader* header= record->GetSimSnarlHeader();

  MCApplication& mcapp = const_cast<MCApplication&>(MCApplication::Instance());
  const VldContext& vldc = *(record -> GetVldContext());
  if ( !(mcapp.GetMC() ) ) {
    // VldContext will be used to construct appropriate geometry only once
    mcapp.InitMC(fMCConfig.c_str(),vldc,fDecayConfig.c_str());
    //ProfilerStart("ptsimprofile");
  }
  mcapp.InitSnarl(vldc);
  
  // Tell the MCApplication Instance about the current User application
  mcapp.SetUserApplication(fMCAppUser);
  
  // Extract StdHep array from record
  // Particle stack will be initialized using StdHep array contents.
  // The particles in this stack represent the initial interaction state 
  // particles + secondaries generated during gminos tracking if a reroot 
  // file is used as input
  TClonesArray* stdheparray = dynamic_cast<TClonesArray*>(const_cast<TObject*>
    (record -> FindComponent("TClonesArray","StdHep")));
 
  if ( !stdheparray || stdheparray->GetEntriesFast() == 0 ) {
    MSG("PTSim",Msg::kWarning) 
    << "No StdHep TClonesArray in SimSnarlRecord or is empty." << endl;
    return JobCResult::kError;
  } 

  // Extract NeuKinList array from record
  // If present, this will be used to determine the number of events
  // in the snarl, which may be more than 1.
  const TClonesArray* neukinarray = dynamic_cast<const TClonesArray*>
    (record -> FindComponent("TClonesArray","NeuKinList"));
  Int_t nevent = 1;
  if ( neukinarray && neukinarray->GetEntriesFast() != 0 ) 
                                nevent = neukinarray->GetEntriesFast();

  // Initialize output DigiScintHit stack, results from tracking will
  // be stored here.
  TClonesArray* hitarray = dynamic_cast<TClonesArray*>(const_cast<TObject*>
    (record -> FindComponent("TClonesArray","DigiScintHits")));
  if ( hitarray ) hitarray -> Clear("C");
  else {
    hitarray = new TClonesArray("DigiScintHit");
    hitarray -> SetName("DigiScintHits");
    record -> AdoptComponent(hitarray);
  }

  // Set random number seeds for this event(snarl) according to the
  // run/snarl number with an optional user specified random seed thrown
  // in as in N. Tagg's DigiPEtoRawDigitModule.
  Int_t run = header -> GetRun();
  Int_t snarl = header -> GetSnarl();
  Int_t seed = run+snarl+fRandomSeed;
  MSG("PTSim",Msg::kDebug) << "Random number seed for run,snarl("
       << run << "," << snarl << ") = " << seed << endl;
  mcapp.GetRandom() -> SetSeed(seed);
  
  fMCAppUser -> SetHitArray(hitarray); // pass the place to store the hits
  // Prime snarl using snarl number, contents of stdheparray, and nevent
  fMCAppUser -> InitSnarl(snarl,stdheparray,nevent); 

  fStopwatch.Start(false);
  // Monte carlo generation takes place here
  gMC -> ProcessRun(nevent);
  fStopwatch.Stop();

  return result;

}

---