MCInitModule.cxx

Go to the documentation of this file.

//
// MCInitModule.cxx
//
// Purpose: A Demo JobControl Module for creating the MC initial state.
//
// Author: S. Kasahara 11/05 
//

#include <iostream>
#include <cassert>
using namespace std;

#include <TClonesArray.h>
#include <TParticle.h>
#include <TSystem.h>
#include <TDatabasePDG.h>

#include "ParticleTransportSim/test/MCInitModule.h"
#include "MessageService/MsgService.h"
#include "JobControl/JobCModuleRegistry.h"
#include "JobControl/JobCommand.h"
#include "Record/SimSnarlRecord.h"
#include "MinosObjectMap/MomNavigator.h"
#include "REROOT_Classes/NeuKin.h"
#include "REROOT_Classes/REROOT_NeuKin.h"

ClassImp(MCInitModule)

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

CVSID("$Id: MCInitModule.cxx,v 1.5 2009/06/22 03:51:01 kasahara Exp $");
JOBMODULE(MCInitModule, "MCInitModule",
                             "A module for initializing the MC snarl.");


// Definition of methods 
// *********************

//......................................................................

const char* MCInitModule::AsString(EInitMethod initmethod) {
  //
  // Purpose: Static method to convert enumerated init method to a string
  //

  switch ( initmethod ) {

  case kSetKine:
    return "SetKine";

  default:
    MSG("MCInit",Msg::kWarning) 
    << "MCInitModule::AsString called with unknown EInitMethod "
    << (int)initmethod << endl;
    return "Unknown";

  } // end of switch
  
}

//......................................................................

MCInitModule::MCInitModule() :
  fRun(0),fSubRun(0),fSnarl(-1),fDetector(Detector::kFar),
  fSimFlag(SimFlag::kMC),fSnarlTimeInterval((time_t)1,900000000),
  fInitMethod(kSetKine),fKineParticle(0) {
  //
  // Purpose: Default constructor
  //
  
}

//......................................................................

MCInitModule::~MCInitModule() {
  // 
  // Purpose:: Destructor
  //

  if ( fKineParticle ) delete fKineParticle; fKineParticle = 0;
  
}

  
//......................................................................

const Registry& MCInitModule::DefaultConfig() const {
  //
  // Purpose: Method to return default configuration.
  // 

  MSG("MCInit",Msg::kDebug) << "MCInitModule::DefaultConfig" << endl;

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

  r.Set("Run", 0);
  r.Set("SubRun", 0);
  r.Set("Snarl", 0);
  r.Set("Detector", (int)Detector::kFar);
  r.Set("SimFlag", (int)SimFlag::kMC);
  r.Set("StartDate",0); 
  r.Set("SnarlTimeInterval",1.9);
  r.Set("InitMethod",(int)kSetKine);
 
  r.LockValues();

  return r;
}

//......................................................................
void MCInitModule::Config(const Registry& r) {
  //
  // Purpose: Configure the module.
  //

  MSG("MCInit",Msg::kDebug) << "MCInitModule::Config" << endl;
  
  Int_t tmpi;
  Double_t tmpd;
  //  const Char_t* tmps;

  if ( r.Get("Run",tmpi) ) fRun = tmpi;
  if ( r.Get("SubRun",tmpi) ) fSubRun = tmpi;
  if ( r.Get("Snarl",tmpi) ) fSnarl = tmpi - 1; // will increment on first Get
  if ( r.Get("SimFlag",tmpi) ) fSimFlag = (SimFlag::SimFlag_t)tmpi;
  if ( r.Get("Detector",tmpi) ) fDetector = (Detector::Detector_t)tmpi;
  
  if ( r.Get("StartDate",tmpi) ) {
    if ( tmpi > 0 ) fStartTimeStamp = VldTimeStamp(tmpi,0,0);
    else fStartTimeStamp = VldTimeStamp(); // set to current time
  }
  if ( r.Get("SnarlTimeInterval",tmpd) ) {
    Int_t sec = (Int_t)tmpd; // truncate
    Int_t nsec = (Int_t)((tmpd - (Double_t)sec)*1.e9);
    fSnarlTimeInterval = VldTimeStamp(sec,nsec);
  }
       
  if ( r.Get("InitMethod",tmpi) ) {
    fInitMethod = (EInitMethod)tmpi;
  }
       
  return;
  
}
    
//......................................................................
void MCInitModule::BeginJob() {
  //
  // Purpose: Called at the beginning of the job, after module configuration
 
  MSG("MCInit",Msg::kDebug) << "MCInitModule::BeginJob " << endl;
  
  // Rewind, because first Get() will move forward
  fStartTimeStamp = fStartTimeStamp - fSnarlTimeInterval; 
  fCurrentVld = VldContext(fDetector,fSimFlag,fStartTimeStamp);

  if ( fInitMethod == kSetKine ) {
    if ( fKineParticle ) {
      MSG("MCInit",Msg::kDebug) << "Initialization method "
          << AsString(kSetKine) << "\npdg " << fKineParticle->GetPdgCode() 
          << "/" << fKineParticle->GetName()  
          << " vtx(x,y,z)(m) (" << fKineParticle->Vx()  << "," 
          <<  fKineParticle->Vy() << "," << fKineParticle->Vz() << ")" 
          << " p(x,y,z)(GeV/c) (" << fKineParticle->Px() << "," 
          <<  fKineParticle->Py() << "," << fKineParticle->Pz() << ")"
          << endl;
    }
    else {
      MSG("MCInit",Msg::kFatal) << "SetKine Initialization method chosen,\n"
          << "but no initial state defined. Use\n"
          << "jc.Mod(''MCInit'').Cmd(''SetKine pdgId vx vy vz px py pz'')\n" 
          << "Abort!" << endl;
      abort();
    }
  }
           
}
  
//......................................................................
JobCResult MCInitModule::Get(MomNavigator *mom) {
  //
  // Purpose:  Create and fill SimSnarlRecord.
  //
  // Notes: A SimSnarlRecord is generated with contents:
  //          i) TClonesArray of name NeuKinList containing objects of type
  //             REROOT_NeuKin.
  //         ii) TClonesArray of name StdHep containing objects of type
  //             TParticle.
  // See description at top of MCInitModule.h for more information.
  //

  JobCResult result(JobCResult::kPassed);  
  MSG("MCInit",Msg::kDebug) << "MCInitModule::Get" << endl;

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

  SimSnarlRecord* simrec = 0;
  // Check for existence of SimSnarlRecord in Mom 
  simrec = dynamic_cast<SimSnarlRecord*>(mom -> GetFragment("SimSnarlRecord"));
  if ( simrec ) {
    MSG("MCInit",Msg::kFatal) << "SimSnarlRecord present in input Mom!" 
                              << endl;
    abort(); // don't want to think about it yet
  }

  fSnarl++;
  fCurrentVld = VldContext(fDetector,fSimFlag,
                           GetCurrentTimeStamp() + fSnarlTimeInterval);
  // Not sure what to do with these flags yet. Copy RerootExodus
  Short_t runtype = 0;
  UInt_t  trigsrc = 0;
  UInt_t errcode = 0;
  // time frames in 1 second buckets, starting with timeframe 0
  Int_t timeframe = (GetCurrentTimeStamp() - fStartTimeStamp 
                   - fSnarlTimeInterval).GetSec(); 
  Int_t spilltype = -1;

  VldTimeStamp now;
  std::string  codename = SimSnarlHeader::TrimCodename("$Name:  $");
  std::string  hostname(gSystem->HostName());

  SimSnarlHeader simheader(fCurrentVld,fRun,fSubRun,runtype,errcode,fSnarl,
                           trigsrc,timeframe,spilltype,
                           now,codename,hostname);
  simrec = new SimSnarlRecord(simheader);

  MSG("MCInit",Msg::kDebug) 
    << "MCInitModule::Get created new SimSnarlRecord.\n" 
    <<  simheader << endl;

  mom -> AdoptFragment(simrec);

  BuildInitialState(simrec);
  
  return result;
  
}

//....................................................................
void MCInitModule::HandleCommand(JobCommand* cmd) {
  // 
  // Handle script commands that don't work through Config:
  //   jc.Mod("MCInit").Cmd("SetKine pdgId vx vy vz px py pz");

  
  if (!cmd->HaveCmd()) return;
  
  string s = cmd->PopCmd();
  if ( s == "SetKine" ) {
    bool isOkay = false;
    if (cmd -> HaveOpt() ) {
      int pdg = cmd->PopIntOpt();
      if ( cmd -> HaveOpt() ) {
        double vx = cmd->PopFloatOpt();
        if ( cmd -> HaveOpt() ) {
          double vy = cmd->PopFloatOpt();
          if ( cmd -> HaveOpt() ) {
            double vz = cmd->PopFloatOpt();
            if ( cmd -> HaveOpt() ) {
              double px = cmd->PopFloatOpt();
              if ( cmd -> HaveOpt() ) {
                double py = cmd -> PopFloatOpt();
                if ( cmd -> HaveOpt() ) {
                  double pz = cmd -> PopFloatOpt();
                  isOkay = true;
                  SetKine(pdg,vx,vy,vz,px,py,pz);
                }
              }
            }
          }
        }
      }
    }
    if (!isOkay ) {
      MSG("MCInit",Msg::kFatal) 
        << "SetKine called w/incorrect parameter list. Abort!" << endl;
      abort();
    }
    return;
  }
    
  // Errors fall through to here
  MSG("MCInit",Msg::kError) << "Illegal command: " << s.c_str() << endl;
  
}

  
//.....................................................................
void MCInitModule::SetKine(Int_t pdg, 
                                Double_t vx, Double_t vy, Double_t vz,
                                Double_t px, Double_t py, Double_t pz) { 
  //
  // Purpose: Method to set initial state kinematics.  
  // 
  // Arguments: particle id (pdg)
  //            vtx & momentum in MINOS standard units
  //            (m and GeV/c respectively)
  // 
  // Notes: fKineParticle (TParticle) will be constructed from the arguments.
  //

  const TDatabasePDG& dbpdg = *(TDatabasePDG::Instance());
  TParticlePDG* partpdg = dbpdg.GetParticle(pdg);

  std::string partname = "???";
  if ( partpdg ) partname=partpdg->GetName();
  
  if ( fKineParticle ) {
    delete fKineParticle; fKineParticle = 0;
  }
  
  // 1 is neugen status code for final state particle 
  // 0 is for neutrino/nucleus primaries
  Int_t status = 1; 
  Int_t mother1 = -1;
  Int_t mother2 = -1;
  Int_t daughter1 = -1;
  Int_t daughter2 = -1;
  Double_t mass = partpdg -> Mass();
  Double_t etot = TMath::Sqrt(px*px + py*py + pz*pz + mass*mass);
  Double_t time = 0;
  fKineParticle = new TParticle(pdg,status,mother1,mother2,daughter1,
                                daughter2,px,py,pz,etot,vx,vy,vz,time);
   
}
  
//......................................................................
void MCInitModule::BuildInitialState(SimSnarlRecord* simrec) {
  //
  //  Purpose: Private method to build initial state particles according
  //           to specified method EInitMethod.  NeuKin and StdHep arrays
  //           are created and filled in SimSnarlRecord.
  //
  //  Arguments: SimSnarlRecord.
  //  
  //  Return: none.
  // 


  switch (fInitMethod) {

  case kSetKine:
    BuildSetKine(simrec);
    break;

  default:
    MSG("MCInit",Msg::kFatal) << "Unknown EInitMethod " 
           << MCInitModule::AsString(fInitMethod) << ". Abort!" << endl;
    abort();
  }

}

//......................................................................
void MCInitModule::BuildSetKine(SimSnarlRecord* simrec) {
  //
  //  Purpose: Private method to build initial state particles according
  //           to specified method EInitMethod::kSetKine.  
  //           NeuKin and StdHep arrays are created and filled in 
  //           SimSnarlRecord.
  //
  //  Arguments: SimSnarlRecord.
  //  
  //  Return: none.
  // 


  // Create StdHep entry
  Int_t nstdhep = 1; // 1 primary track
  TClonesArray* stdheparray = new TClonesArray("TParticle",nstdhep);
  stdheparray -> SetName("StdHep");
  simrec -> AdoptComponent(stdheparray);

  if ( !fKineParticle ) {
    MSG("MCInit",Msg::kFatal) << "SetKine Initialization method chosen,\n"
          << "but no initial state defined. Use\n"
          << "jc.Mod(''MCInit'').Cmd(''SetKine pdgId vx vy vz px py pz'')\n" 
          << "Abort!" << endl;
    abort();
  }
  new ((*stdheparray)[0])TParticle(*fKineParticle);


  // Create single dummy REROOT_NeuKin object, REROOT_NeuKin to be replaced
  Int_t nneukin = 1; // 1 primary
  TClonesArray* neukinarray = new TClonesArray("REROOT_NeuKin",nneukin);
  neukinarray -> SetName("NeuKinList");
  simrec -> AdoptComponent(neukinarray);

  // This bit of nonsense is because REROOT_NeuKin constructor doesn't
  // initialize data members
  NEUKIN_DEF neukin;
  neukin.ID= 1; //index of this neukin, set to 1 cuz REROOT_NeuKin subtracts 1
  neukin.INu = 0; // particle id
  neukin.INuNoOsc = 0;
  neukin.ITg = 0;
  neukin.IBoson = -99999; // to be the same as reroot
  neukin.IResonance = 0;
  neukin.IAction = 0;
  neukin.A = 0;
  neukin.Z = 0;
  neukin.Sigma = 0;
  for ( int ic = 0; ic < 4; ic++ ) {
    neukin.P4Neu[ic] = 0;
    neukin.P4NeuNoOsc[ic] = 0;
    neukin.P4Tgt[ic] = 0;
  }
  neukin.X = 0;
  neukin.Y = 0;
  neukin.Q2 = 0;
  neukin.W2 = 0;
  for ( int ic = 0; ic < 4; ic++ ) {
    neukin.P4Shw[ic] = 0;
    neukin.P4Mu1[ic] = 0;
    neukin.P4Mu2[ic] = 0;
  }
  if (TMath::Abs(fKineParticle->GetPdgCode()) == 13 ) {
    // Muon primary, gminos/reroot labels these muon final state
    neukin.P4Mu1[0] = fKineParticle->Px();
    neukin.P4Mu1[1] = fKineParticle->Py();
    neukin.P4Mu1[2] = fKineParticle->Pz();
    neukin.P4Mu1[3] = fKineParticle->Energy();
  }
  neukin.EMFrac = 0;
  for ( int ic = 0; ic < 4; ic++ ) {
    neukin.P4El1[ic] = 0;
    neukin.P4El2[ic] = 0;
    neukin.P4Tau[ic] = 0;
  }
  
  // new with placement, create dummy REROOT_NeuKin
  new ((*neukinarray)[0])REROOT_NeuKin(&neukin);

  
  return;
  
}

---