LightInjector.cxx

Go to the documentation of this file.

// $Id: LightInjector.cxx,v 1.9 2006/12/12 17:07:04 rhatcher Exp $
//
// Job control interface to input data streams
//
// messier@huhepl.harvard.edu
#include "TSystem.h"
#include "TRegexp.h"

#include "LightInjector.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 "Digitization/DigiList.h"
#include "Digitization/DigiPE.h"
#include "Digitization/DigiSignal.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 "Conventions/Munits.h"

#include <algorithm>

CVSID("$Id: LightInjector.cxx,v 1.9 2006/12/12 17:07:04 rhatcher Exp $");
JOBMODULE(LightInjector,"LightInjector","Simulated Light Injection");
PlexPixelSpotId LightInjector::fPsid;
TFile*          LightInjector::fFile=0;
TH1*            LightInjector::fHisto[65]; // One for each pixel.


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

LightInjector::LightInjector() : 
  fCurrentRun(1000),
  fCurrentSnarl(0)
{ }

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

LightInjector::~LightInjector() 
{ 
}

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

void LightInjector::BeginJob() 
{
}


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

const Registry& LightInjector::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("plane",        2);
  r.Set("strip",        10);
  r.Set("end",          2);
  r.Set("pulseHeight",  10.0);
  r.Set("pulseWidth",   0.0);
  r.LockValues();

  return r;
}

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

void LightInjector::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; }
  if (r.Get("plane",tmpi)) { fPlane = tmpi; }
  if (r.Get("strip",tmpi)) { fStrip = tmpi; }
  if (r.Get("end",tmpi)) { fEnd = tmpi; }
  if (r.Get("pulseHeight",tmpd)) { fPulseHeight = tmpd; }
  if (r.Get("pulseWidth",tmpd)) { fPulseWidth = tmpd; }

  fPulseTime = VldTimeStamp(fDate,fTime,0),
  fContext =   VldContext(Detector::CharToEnum((char)fDetector),
                          SimFlag::kReroot,
                          fPulseTime),
  fPulseSeid = PlexStripEndId(Detector::CharToEnum((char)fDetector),
                          fPlane,
                          fStrip,
                          StripEnd::CharToEnum((char)fEnd));
}

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

JobCResult LightInjector::Get()
{
  MomNavigator* mom = this->GetMom();
  assert(mom);
  mom -> Clear(); // Moving on so clear contents of Mom


  // Make a SimSnarl.
  Int_t   trigbits = 0;
  Short_t subrun = 0;
  Short_t runtype = 0;
  Int_t   errcode = 0;
  Int_t   timeframe = -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);
  fPulseTime = VldTimeStamp(fDate,fTime,0);
  fPulseTime.Add(offset);

  // Reset the context.
  fContext = VldContext(Detector::CharToEnum((char)fDetector),
                        SimFlag::kMC,
                        fPulseTime);

  // Get a PlexHandle to find the psid.
  PlexHandle plex(fContext);
  fPsid = plex.GetPixelSpotId(fPulseSeid);
  if(!fPsid.IsValid()) {
    MSG("DetSim",Msg::kError) << "Invalid psid " << fPsid.AsString()
                              << " for strip " << fPulseSeid.AsString()
                              << " at context " << fContext.AsString()
                              << endl;
    return JobCResult::kError;
  }

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

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

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

  // Add the DigiPE list to the simsnarl.
  TObjArray* PeList = new TObjArray(0);
  PeList->SetName("DigiListPe");
  PeList->SetOwner(true);

  static double tot_pulses = 0;
  static double tot_npe = 0;
  
  // Find out how many of them to make.
  gRandom->SetSeed(fCurrentSnarl);
  Int_t npe = gRandom->Poisson(fPulseHeight);

  // Make the photoelectron objects.
  for(int i=0;i<npe;i++) {
    // Time is 100ns +- pulse width.
    double time = gRandom->Gaus(100.0*Munits::ns, fPulseWidth);
    
    DigiPE*  pe = new DigiPE(time, fPsid, DigiSignal::kGenuine);
    PeList->Add(pe);
  }

  tot_pulses +=1;
  tot_npe += npe;
  MSG("DetSim",Msg::kInfo) << "SnarlNo: " << fCurrentSnarl << "  Pulse: " << npe << "avg: " << tot_npe/tot_pulses << " pe into " << fPsid.AsString()  << endl;

  simsnarl->AdoptTemporary(PeList);

  return JobCResult::kAOK;
}

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

JobCResult LightInjector::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 LightInjector::Prev(int n) 
{
  if(fCurrentSnarl>0) fCurrentSnarl-=n;
  return this->Get();
}

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

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

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

JobCResult LightInjector::Ana(const MomNavigator* mom)
{
  // Makes histograms of all the pixels
  // in the pulsed PMT.
  MSG("DetSim",Msg::kVerbose) << "LI Ana:" << endl;
  
  if(fFile==0) {
    TFile* save = gFile;
    fFile = new TFile("mcli.root","RECREATE");
    fFile->cd();
    for(int c=0;c<65;c++) {
      fHisto[c] = new TH1F(Form("h%02d",c),
                           Form("Pixel %02d",c),
                           20500,-500,20000 );
    }
    if(save) save->cd();
  }
 

  int ndigits_tot = 0;
  int ndigits = 0;

  // Get the RawDigitDataBlock.
  TObject* tobj;
  TIter    fragiter = mom->FragmentIter();
  while( ( tobj = fragiter.Next() ) ) {
    RawRecord* rawRec = dynamic_cast<RawRecord*>(tobj);
    if(rawRec) {
      TIter itr = rawRec->GetRawBlockIter();
      RawDataBlock* rawBlk;
      while ( ( rawBlk = dynamic_cast<RawDataBlock*>(itr.Next()) ) ) {
        RawDigitDataBlock* rddb = dynamic_cast<RawDigitDataBlock*>(rawBlk);
        if(rddb) {
          
          cout << "Ana context: " << rawRec->GetRawHeader()->GetVldContext().AsString() << endl;
          PlexHandle plex(rawRec->GetRawHeader()->GetVldContext());

          int npixels = 16;
          if(fPsid.GetElecType()==ElecType::kQIE) npixels = 64;
          
          PlexPixelSpotId pixels[65];
          RawChannelId    channels[65];
          Double_t        values[65];
          for(int i=0;i<npixels;i++){
            pixels[i] = fPsid;
            pixels[i].SetPixel(i);
            channels[i] = plex.GetRawChannelId(pixels[i]);
            values[i]=0;
          }

          TIter iter = rddb->GetDatumIter();
          while( (tobj=iter.Next() ) ) {
            
            // Process a digit.     
            RawDigit*        rawDigit = dynamic_cast<RawDigit *>(tobj);
            
            if(rawDigit) {
              // Try to add it.
              // We do it this way to 
              // make sure we get all buckets.
              ndigits_tot++;
              for(int i=0;i<npixels;i++) {
                if(rawDigit->GetChannel().IsSameChannel(channels[i])) {
                  ndigits++;
                  values[i] += rawDigit->GetADC();
                  if(rawDigit->GetChannel().GetElecType()==ElecType::kQIE)
                    values[i] -= 50; // Remove digital offset for each bucket.
                }
              }    
            }

          }


          MSG("DetSim",Msg::kDebug) << "Got " << ndigits << " digits; total in event = " << ndigits_tot << endl;
          // Fill histograms.
          for(int i=0;i<npixels;i++) {
            fHisto[i]->Fill(values[i]);
          }
          


        }       
      }
    }
  } // It's so ugly!
  

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

---