PhotonInjector.cxx

Go to the documentation of this file.

// $Id: PhotonInjector.cxx,v 1.6 2007/11/11 05:47:55 rhatcher Exp $
//
// FILL_IN: [Document your code!!]
//
// n.tagg1@physics.ox.ac.uk
#include "PhotonInjector.h"

#include <TRandom.h>
#include <TRandom3.h>

#include "MessageService/MsgService.h"
#include "MinosObjectMap/MomNavigator.h"
#include "JobControl/JobCModuleRegistry.h" // For JOBMODULE macro

#include "Digitization/DigiScintHit.h"
#include "Digitization/DigiPE.h"
#include "Digitization/DigiSignal.h"

#include "Record/SimSnarlRecord.h"
#include "Record/SimSnarlHeader.h"

#include "PhotonTransportModule.h"
#include "PhotonTransportMaker.h"
#include "PhotonConfiguration.h"

#include "Plex/PlexHandle.h"

#include "Calibrator/Calibrator.h"

#include <TFile.h>
#include <TNtuple.h>
#include <TSystem.h>

#include <cmath>

JOBMODULE(PhotonInjector, "PhotonInjector",
          "Photon transport from light injection to photoelectroncs");
CVSID("$Id: PhotonInjector.cxx,v 1.6 2007/11/11 05:47:55 rhatcher Exp $");


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

PhotonInjector::PhotonInjector() :
  fRandom(new TRandom3),
  fRandomSeed(0)
{
}

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

PhotonInjector::~PhotonInjector() 
{
}

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

JobCResult PhotonInjector::Get()
{
  //======================================================================
  // Does the work. 
  // Creates photons near the LIM ashtray
  //======================================================================
  MomNavigator* mom = this->GetMom();
  assert(mom);
  mom -> Clear(); // Moving on so clear contents of Mom
  
  MSG("Photon",Msg::kInfo) << "Simulated Light Injection  running." << endl;
  
  // Make a SimSnarl.
  Int_t   trigbits = 0x2; // RawLI
  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(fCurrentSnarl/fRate,&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,fRun,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);


  TIter    fragiter = mom->FragmentIter();

  // Set the random seed to match the event.
  // Also add some user-configurable randomness, in case they want to
  // generate the same event with a different outcome.
  fRandom->SetSeed( fRun 
                    + fCurrentSnarl
                    + fRandomSeed );
  
  // We need a context to do table lookups
  VldContext simContext = fContext;

  // But.. we want our context to be "MC" not "Reroot". Ensure that's correct.
  fContext = VldContext(simContext.GetDetector(),SimFlag::kMC,simContext.GetTimeStamp());

  Calibrator::Instance().Reset(fContext);

  // Create the stat report.
  JobCResult res;

  fPeList = new TObjArray(0);
  fPeList->SetName("DigiListPe");
  fPeList->SetOwner(true);
  
  //
  // Do the actual work.
  //
  SimulateEvent();

  // Store the results.
  simsnarl->AdoptTemporary(fPeList);
  
  
  return JobCResult::kPassed; // kNoDecision, kFailed, etc.
}

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

const Registry& PhotonInjector::DefaultConfig() const
{
//======================================================================
// Supply the default configuration for the module
//======================================================================
  static Registry kReg;
  kReg = PhotonConfiguration();
  kReg.Set("Date",VldTimeStamp().GetDate()); // Now
  kReg.Set("Time",VldTimeStamp().GetTime()); // Now
  kReg.Set("Detector",2); // Far
  kReg.Set("PulserBox",1);
  kReg.Set("Led",1);
  kReg.Set("Rate",200.); //Hz
  kReg.Set("PulseWidth",50e-9);
  kReg.Set("PulseHeight",330.);
  return kReg;
}

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

void PhotonInjector::Config(const Registry& r)
{
//======================================================================
// Configure the module given the Registry r
//======================================================================
  r.Get("Run",fRun);
  r.Get("Date",fDate);
  r.Get("Time",fTime);
  r.Get("Detector",fDetector);
  r.Get("PulserBox",fPulserBox);
  r.Get("Led",fLed);
  r.Get("Rate",fRate);
  r.Get("PulseWidth",fPulseWidth);
  r.Get("PulseHeight",fPulseHeight);
  r.Get("RandomSeed",fRandomSeed);
}


JobCResult PhotonInjector::SimulateEvent()
{
  JobCResult res = JobCResult::kPassed;


  PlexHandle plex(fContext);
  UgliGeomHandle ugli(fContext);

  // Identify the LED.
  PlexLedId led(fContext.GetDetector(), fPulserBox, fLed);
  std::vector<PlexStripEndId> stripsHit = plex.GetStripEndIdVector(led);

  cout << "Pulsing LED: " << led.AsString() << endl;
  cout << "Total stripends = " << stripsHit.size() << endl;

  // Loop over the strip ends.
  for(UInt_t i=0; i<stripsHit.size(); i++) {
    PlexStripEndId seid = stripsHit[i];
    PlexStripEndId farseid = seid;
    farseid.SetEnd((seid.GetEnd()==StripEnd::kNegative)?StripEnd::kPositive:StripEnd::kNegative);
    
    // Total length of green fibre
    UgliStripHandle ustrip = ugli.GetStripHandle(seid);
    Float_t xfibre = ustrip.GetHalfLength() + ustrip.WlsPigtail(seid.GetEnd());
    if(seid.GetEnd()==StripEnd::kNegative) xfibre = -xfibre;

    // Clear fibre is the same on both ends, and can be ignored.

    // Find attenuation factor.
    Float_t atten = Calibrator::Instance().DecalAttenCorrected(1,xfibre,seid);

    // Find number of pe at each end.
    Int_t npe_near = fRandom->Poisson(fPulseHeight);
    Int_t npe_far  = fRandom->Poisson(fPulseHeight/atten);
    PlexPixelSpotId spot_near = plex.GetPixelSpotId(seid);
    PlexPixelSpotId spot_far  = plex.GetPixelSpotId(farseid);
    std::cout << "Generating near/far: " << npe_near << "/" << npe_far << std::endl;

    for(int i=0;i<(npe_near+npe_far);i++) {
      PlexPixelSpotId spot = spot_near;
      if(i>=npe_near) spot = spot_far;
      if(!spot.IsValid()) continue;
      
      // Choose time.
      Double_t t = gRandom->Gaus()*fPulseWidth
        + fContext.GetTimeStamp().GetNanoSec()*1e-9;
      
      DigiPE* pe = new DigiPE(t,spot,DigiSignal::kFibreLight); // Close enough.
      fPeList->Add(pe);
    }
     
  }

  return res;
}

---