ToyMCModule.cxx

Go to the documentation of this file.

// $Id: ToyMCModule.cxx,v 1.14 2007/11/11 06:16:02 rhatcher Exp $
//

#include "TSystem.h"
#include "TRegexp.h"

#include "ToyMCModule.h"
#include <cassert>
#include <cmath>
#include "MessageService/MsgService.h"
#include "MinosObjectMap/MomNavigator.h"
#include "JobControl/JobCInputModule.h"
#include "JobControl/JobCModuleRegistry.h"
#include "JobControl/JobCEnv.h"
#include "RawData/RawRecord.h"
#include "RawData/RawDaqSnarlHeader.h"
#include "Record/RecRecord.h"
#include "Record/RecPhysicsHeader.h"
#include "Registry/Registry.h"
#include "Validity/VldContext.h"
#include "Util/UtilString.h"
#include <TRandom.h>
#include "UgliGeometry/UgliGeomHandle.h"
#include "Digitization/DigiList.h"
#include "Digitization/DigiScintHit.h"
#include "Record/SimSnarlRecord.h"
#include "Record/SimSnarlHeader.h"
#include "RawData/RawRecord.h"
#include "RawData/RawDaqSnarlHeader.h"
#include "RawData/RawDigitDataBlock.h"
#include "RawData/RawDigit.h"
#include <TClonesArray.h>
#include <TParticle.h>
#include "Conventions/Munits.h"
#include "Util/LoadMinosPDG.h"

#include "DataUtil/Truthifier.h"

#include <algorithm>

CVSID("$Id: ToyMCModule.cxx,v 1.14 2007/11/11 06:16:02 rhatcher Exp $");
JOBMODULE(ToyMCModule,"ToyMCModule","Toy Muon MC");


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

ToyMCModule::ToyMCModule() : 
  fCurrentSnarl(0)
{
  LoadMinosPDG();
}

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

ToyMCModule::~ToyMCModule() 
{ 
}

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

void ToyMCModule::BeginJob() 
{
}


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

const Registry& ToyMCModule::DefaultConfig() const 
{
//======================================================================
// Get the default configuration for this module
//======================================================================
  static Registry r; // Default configuration for module

  // Set name of config
  std::string name = this->GetName();
  name += ".config.default";
  r.SetName(name.c_str());

  // Set values in configuration
  r.UnLockValues();
  r.Set("DetectorType", 2);
  r.Set("Rate",         1.0);
  r.Set("Date",         20030101);
  r.Set("Time",         000000);
  r.Set("StartPlane",   1);
  r.Set("StopPlane",    50);
  r.Set("StartX",       0.30*Munits::m);
  r.Set("StartY",       0.30*Munits::m);
  r.Set("DX",       0.*Munits::m);
  r.Set("DY",       0.*Munits::m);

  r.Set("EnergyPerPlane",2.0*Munits::MeV);
  r.LockValues();

  return r;
}

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

void ToyMCModule::Config(const Registry& r) 
{
//======================================================================
// Configure the module given the Registry r
//======================================================================
  int    tmpi;
  double tmpd;

  if (r.Get("DetectorType",tmpi)) { fDetector = tmpi; }
  if (r.Get("Rate",tmpd)) { fRate = tmpd; }
  if (r.Get("Date",tmpi)) { fDate = tmpi; }
  if (r.Get("Time",tmpi)) { fTime = tmpi; }
  
  r.Get("StartPlane",fStartPlane);
  r.Get("StopPlane",fStopPlane);
  r.Get("StartX",fStartX);
  r.Get("StartY",fStartY);
  r.Get("DX",fDX);
  r.Get("DY",fDY);
  r.Get("EnergyPerPlane",fEnergyPerPlane);
}

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

JobCResult ToyMCModule::Get()
{
  MomNavigator* mom = this->GetMom();
  assert(mom);
  mom -> Clear(); // Moving on so clear contents of Mom
  
  MSG("ToyMC",Msg::kInfo) << "Toy MC running." << endl;
  
  // Make a SimSnarl.
  Int_t   run      = 1000;
  Int_t   trigbits = 0;
  Short_t subrun = 0;
  Short_t runtype = 0;
  Int_t   errcode = 0;
  Int_t   timeframe = -1;
  Int_t   spilltype = -1;
  
  // Increment the time by the pulse rate.
  int secs;
  double frac = frexp(fRate*fCurrentSnarl,&secs);
  int nsecs = int(frac*1e9);
  VldTimeStamp offset(secs,nsecs);
  VldTimeStamp fPulseTime = VldTimeStamp(fDate,fTime,0);
  fPulseTime.Add(offset);
  
  // Reset the context.
  fContext = VldContext(Detector::CharToEnum((char)fDetector),
                        SimFlag::kMC,
                        fPulseTime);
  
  VldTimeStamp now;
  std::string  codename = SimSnarlHeader::TrimCodename("$Name:  $");
  std::string  hostname(gSystem->HostName());

  SimSnarlHeader simheader(fContext,run,subrun,runtype,
                           errcode,fCurrentSnarl,trigbits,timeframe,spilltype,
                           now,codename,hostname);
                                                
  SimSnarlRecord* simsnarl = new SimSnarlRecord(simheader);

  simsnarl->GetTempTags().Set("stream","SimSnarl");  // no idea....
    
  // Give it to mom.
  mom -> AdoptFragment(simsnarl);

  // Generate the event:
  gRandom->SetSeed(fCurrentSnarl);


  UgliGeomHandle ugli(fContext);  

  std::vector<DigiScintHit> hits;

  for(int plane=fStartPlane;plane<=fStopPlane;plane++) {
    PlexPlaneId planeId(Detector::CharToEnum((char)fDetector),plane);
    if(!(planeId.IsValid())) break;
    
    UgliScintPlnHandle uplane = ugli.GetScintPlnHandle(planeId);
    if(!(uplane.IsValid())) break;
    std::vector<UgliStripHandle> strips = uplane.GetStripHandleVector();
    
    
    // Find the strip that best matches our coordinates.
    double min = 1e19;
    UgliStripHandle ustrip;

    TVector3 x(fStartX + (double)(fCurrentSnarl)*fDX,
               fStartY + (double)(fCurrentSnarl)*fDY,
               uplane.GetCenter().z());


    for(UInt_t istrip = 0;istrip<strips.size(); istrip++) {
      if(strips[istrip].IsValid()) {
        double dist = fabs(strips[istrip].GlobalToLocal(x).y());
        if(dist<min) { 
          min = dist; 
          ustrip = strips[istrip];
        }
      }
    }
    if(!ustrip.IsValid()) break;

    // Find the coords on the strip.
    TVector3 x0 = ustrip.GlobalToLocal(x);

    MSG("ToyMC",Msg::kInfo) << "Plane " << plane
                            << "  x: (" << x.x() <<","<<  x.y()  <<","<<  x.z() << ")\t"
                            << "x0: (" << x0.x() <<","<<  x0.y()  <<","<<  x0.z() << ")\t"
      //<< "x1: (" << x1.x() <<","<<  x1.y()  <<","<<  x1.z() << ")\t"
      //<< "x2: (" << x2.x() <<","<<  x2.y()  <<","<<  x2.z() << ")"
                            << endl;

    // Ensure we're in the scintillator.
    if( x0.y() > ustrip.GetHalfWidth()) {
      x0.SetY(ustrip.GetHalfWidth()-1.0*Munits::mm);
    }

    if( x0.y() < -ustrip.GetHalfWidth()) {
      x0.SetY(-ustrip.GetHalfWidth()+1.0*Munits::mm);
    }


    if( fabs(x0.x()) > ustrip.GetHalfLength() ) {
      MSG("ToyMC",Msg::kWarning) << "Outside strip length!" << endl;
      break;
    }
    TVector3 x1 = x0; 
    x1.SetZ(-(ustrip.GetHalfThickness()));
    TVector3 x2 = x0; 
    x2.SetZ( (ustrip.GetHalfThickness()));

    double t1 = (x.z() - (ustrip.GetHalfThickness()))/Munits::c_light;
    double t2 = (x.z() + (ustrip.GetHalfThickness()))/Munits::c_light;

    DigiScintHit hit(1,  // track ID <<","<< 
                     13, // muon,
                     1e9,
                     ustrip.GetSEId(),
                     t1,
                     x1.x(), x1.y(), x1.z(),
                     t2,
                     x2.x(), x2.y(), x2.z(),
                     (ustrip.GetHalfThickness())*2,
                     fEnergyPerPlane);
    MSG("ToyMC",Msg::kVerbose) << hit;
    hits.push_back(hit);
                                         
  }

  TClonesArray *hitlist = new TClonesArray("DigiScintHit",hits.size());
  hitlist->SetName("DigiScintHits");
  for(unsigned int i=0;i<hits.size();i++) {
    new ((*hitlist)[i])  DigiScintHit( hits[i] );
  }
  simsnarl->AdoptComponent(hitlist);

  // Add truth.
  // New TClonesArray of the right size
  TClonesArray* newstdhep = new TClonesArray("TParticle",2);
  newstdhep->SetName("StdHep");
  new((*newstdhep)[0])
    TParticle(13, // muon
              0, 0, 0, 0, 0,
              0.,0.,10000.,10000.,
              0.,0.,1.0,0);
  new((*newstdhep)[1])
    TParticle(13, // muon
              0, 0, 0, 0, 0,
              0.,0.,10000.,10000.,
              0.,0.,1.0,0);
  
  simsnarl->AdoptComponent(newstdhep);

  return JobCResult::kAOK;
}

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

JobCResult ToyMCModule::Next(int n)
{
//======================================================================
// Advance the position in the stream n record sets. Load the records
// at the last position
//======================================================================
  fCurrentSnarl+=n;
  return this->Get();
}

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

JobCResult ToyMCModule::Prev(int n) 
{
  if(fCurrentSnarl>0) fCurrentSnarl-=n;
  return this->Get();
}

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

JobCResult ToyMCModule::GoTo(int /* run */, int snarl, int /* searchDir */) 
{
  fCurrentSnarl = snarl;
  return this->Get();
}

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

JobCResult ToyMCModule::Ana(const MomNavigator* mom)
{
  //const Truthifier& truth = Truthifier::Instance(mom);

  // Find PrimaryCandidateRecord fragment in MOM.
  CandRecord *candrec = dynamic_cast<CandRecord *>
    (mom->GetFragment("CandRecord", "PrimaryCandidateRecord"));
  if (candrec == 0) {
    MSG("toymc", Msg::kWarning) << "No PrimaryCandidateRecord in  MOM."
                                  << endl;
    JobCResult result;
    result.SetWarning().SetFailed();
    return result;
  }

  // Get CandDigits.
  CandDigitListHandle *cdlh_ptr = dynamic_cast<CandDigitListHandle *>
    (candrec->FindCandHandle("CandDigitListHandle", "canddigitlist"));
  
   if(cdlh_ptr==0) return JobCResult::kError;
  CandDigitListHandle cdlh = *cdlh_ptr;

  CandDigitHandleItr cdhItr(cdlh.GetDaughterIterator());
  //while(CandDigitHandle *cdh = cdhItr()) {
    //cout << truth.BestSEIdOfDigit(*cdh).AsString() << endl;
  //}


  return JobCResult::kPassed; // kNoDecision, kFailed, etc.
}

---