loon: LightInjector.cxx Source File

00001 
00002 // $Id: LightInjector.cxx,v 1.9 2006/12/12 17:07:04 rhatcher Exp $
00003 //
00004 // Job control interface to input data streams
00005 //
00006 // messier@huhepl.harvard.edu
00008 #include "TSystem.h"
00009 #include "TRegexp.h"
00010 
00011 #include "LightInjector.h"
00012 #include <cassert>
00013 #include <cmath>
00014 
00015 #include "MessageService/MsgService.h"
00016 #include "MinosObjectMap/MomNavigator.h"
00017 #include "JobControl/JobCInputModule.h"
00018 #include "JobControl/JobCModuleRegistry.h"
00019 #include "JobControl/JobCEnv.h"
00020 #include "RawData/RawRecord.h"
00021 #include "RawData/RawDaqSnarlHeader.h"
00022 #include "Record/RecRecord.h"
00023 #include "Record/RecPhysicsHeader.h"
00024 #include "Registry/Registry.h"
00025 #include "Validity/VldContext.h"
00026 #include "Util/UtilString.h"
00027 #include <TRandom.h>
00028 #include "Digitization/DigiList.h"
00029 #include "Digitization/DigiPE.h"
00030 #include "Digitization/DigiSignal.h"
00031 #include "Record/SimSnarlRecord.h"
00032 #include "Record/SimSnarlHeader.h"
00033 #include "RawData/RawRecord.h"
00034 #include "RawData/RawDaqSnarlHeader.h"
00035 #include "RawData/RawDigitDataBlock.h"
00036 #include "RawData/RawDigit.h"
00037 
00038 #include "Conventions/Munits.h"
00039 
00040 #include <algorithm>
00041 
00042 CVSID("$Id: LightInjector.cxx,v 1.9 2006/12/12 17:07:04 rhatcher Exp $");
00043 JOBMODULE(LightInjector,"LightInjector","Simulated Light Injection");
00044 PlexPixelSpotId LightInjector::fPsid;
00045 TFile*          LightInjector::fFile=0;
00046 TH1*            LightInjector::fHisto[65]; // One for each pixel.
00047 
00048 
00049 //......................................................................
00050 
00051 LightInjector::LightInjector() : 
00052   fCurrentRun(1000),
00053   fCurrentSnarl(0)
00054 { }
00055 
00056 //......................................................................
00057 
00058 LightInjector::~LightInjector() 
00059 { 
00060 }
00061 
00062 //......................................................................
00063 
00064 void LightInjector::BeginJob() 
00065 {
00066 }
00067 
00068 
00069 //......................................................................
00070 
00071 const Registry& LightInjector::DefaultConfig() const 
00072 {
00073 //======================================================================
00074 // Get the default configuration for this module
00075 //======================================================================
00076   static Registry r; // Default configuration for module
00077 
00078   // Set name of config
00079   std::string name = this->GetName();
00080   name += ".config.default";
00081   r.SetName(name.c_str());
00082 
00083   // Set values in configuration
00084   r.UnLockValues();
00085   r.Set("detectorType", 2);
00086   r.Set("rate",         1.0);
00087   r.Set("date",         20030101);
00088   r.Set("time",         000000);
00089   r.Set("plane",        2);
00090   r.Set("strip",        10);
00091   r.Set("end",          2);
00092   r.Set("pulseHeight",  10.0);
00093   r.Set("pulseWidth",   0.0);
00094   r.LockValues();
00095 
00096   return r;
00097 }
00098 
00099 //......................................................................
00100 
00101 void LightInjector::Config(const Registry& r) 
00102 {
00103 //======================================================================
00104 // Configure the module given the Registry r
00105 //======================================================================
00106   int    tmpi;
00107   double tmpd;
00108 
00109   if (r.Get("detectorType",tmpi)) { fDetector = tmpi; }
00110   if (r.Get("rate",tmpd)) { fRate = tmpd; }
00111   if (r.Get("date",tmpi)) { fDate = tmpi; }
00112   if (r.Get("time",tmpi)) { fTime = tmpi; }
00113   if (r.Get("plane",tmpi)) { fPlane = tmpi; }
00114   if (r.Get("strip",tmpi)) { fStrip = tmpi; }
00115   if (r.Get("end",tmpi)) { fEnd = tmpi; }
00116   if (r.Get("pulseHeight",tmpd)) { fPulseHeight = tmpd; }
00117   if (r.Get("pulseWidth",tmpd)) { fPulseWidth = tmpd; }
00118 
00119   fPulseTime = VldTimeStamp(fDate,fTime,0),
00120   fContext =   VldContext(Detector::CharToEnum((char)fDetector),
00121                           SimFlag::kReroot,
00122                           fPulseTime),
00123   fPulseSeid = PlexStripEndId(Detector::CharToEnum((char)fDetector),
00124                           fPlane,
00125                           fStrip,
00126                           StripEnd::CharToEnum((char)fEnd));
00127 }
00128 
00129 //......................................................................  
00130 
00131 JobCResult LightInjector::Get()
00132 {
00133   MomNavigator* mom = this->GetMom();
00134   assert(mom);
00135   mom -> Clear(); // Moving on so clear contents of Mom
00136 
00137 
00138   // Make a SimSnarl.
00139   Int_t   trigbits = 0;
00140   Short_t subrun = 0;
00141   Short_t runtype = 0;
00142   Int_t   errcode = 0;
00143   Int_t   timeframe = -1;
00144 
00145   // Increment the time by the pulse rate.
00146   int secs;
00147   double frac = frexp(fRate*fCurrentSnarl,&secs);
00148   int nsecs = int(frac*1e9);
00149   VldTimeStamp offset(secs,nsecs);
00150   fPulseTime = VldTimeStamp(fDate,fTime,0);
00151   fPulseTime.Add(offset);
00152 
00153   // Reset the context.
00154   fContext = VldContext(Detector::CharToEnum((char)fDetector),
00155                         SimFlag::kMC,
00156                         fPulseTime);
00157 
00158   // Get a PlexHandle to find the psid.
00159   PlexHandle plex(fContext);
00160   fPsid = plex.GetPixelSpotId(fPulseSeid);
00161   if(!fPsid.IsValid()) {
00162     MSG("DetSim",Msg::kError) << "Invalid psid " << fPsid.AsString()
00163                               << " for strip " << fPulseSeid.AsString()
00164                               << " at context " << fContext.AsString()
00165                               << endl;
00166     return JobCResult::kError;
00167   }
00168 
00169   VldTimeStamp now;
00170   std::string  codename = SimSnarlHeader::TrimCodename("$Name:  $");
00171   std::string  hostname(gSystem->HostName());
00172 
00173   SimSnarlHeader simheader(fContext,fCurrentRun,subrun,runtype,
00174                            errcode,fCurrentSnarl,trigbits,timeframe,-1,
00175                            now,codename,hostname);
00176                                                 
00177   SimSnarlRecord* simsnarl = new SimSnarlRecord(simheader);
00178 
00179   simsnarl->GetTempTags().Set("stream","SimSnarl");  // no idea....
00180     
00181   // Give it to mom.
00182   mom -> AdoptFragment(simsnarl);
00183 
00184   // Add the DigiPE list to the simsnarl.
00185   TObjArray* PeList = new TObjArray(0);
00186   PeList->SetName("DigiListPe");
00187   PeList->SetOwner(true);
00188 
00189   static double tot_pulses = 0;
00190   static double tot_npe = 0;
00191   
00192   // Find out how many of them to make.
00193   gRandom->SetSeed(fCurrentSnarl);
00194   Int_t npe = gRandom->Poisson(fPulseHeight);
00195 
00196   // Make the photoelectron objects.
00197   for(int i=0;i<npe;i++) {
00198     // Time is 100ns +- pulse width.
00199     double time = gRandom->Gaus(100.0*Munits::ns, fPulseWidth);
00200     
00201     DigiPE*  pe = new DigiPE(time, fPsid, DigiSignal::kGenuine);
00202     PeList->Add(pe);
00203   }
00204 
00205   tot_pulses +=1;
00206   tot_npe += npe;
00207   MSG("DetSim",Msg::kInfo) << "SnarlNo: " << fCurrentSnarl << "  Pulse: " << npe << "avg: " << tot_npe/tot_pulses << " pe into " << fPsid.AsString()  << endl;
00208 
00209   simsnarl->AdoptTemporary(PeList);
00210 
00211   return JobCResult::kAOK;
00212 }
00213 
00214 //......................................................................
00215 
00216 JobCResult LightInjector::Next(int n)
00217 {
00218 //======================================================================
00219 // Advance the position in the stream n record sets. Load the records
00220 // at the last position
00221 //======================================================================
00222   fCurrentSnarl+=n;
00223   return this->Get();
00224 }
00225 
00226 //......................................................................
00227 
00228 JobCResult LightInjector::Prev(int n) 
00229 {
00230   if(fCurrentSnarl>0) fCurrentSnarl-=n;
00231   return this->Get();
00232 }
00233 
00234 //......................................................................
00235 
00236 JobCResult LightInjector::GoTo(int /* run */, int snarl, int /* searchDir */) 
00237 {
00238   fCurrentSnarl = snarl;
00239   return this->Get();
00240 }
00241 
00242 //......................................................................
00243 
00244 JobCResult LightInjector::Ana(const MomNavigator* mom)
00245 {
00246   // Makes histograms of all the pixels
00247   // in the pulsed PMT.
00248   MSG("DetSim",Msg::kVerbose) << "LI Ana:" << endl;
00249   
00250   if(fFile==0) {
00251     TFile* save = gFile;
00252     fFile = new TFile("mcli.root","RECREATE");
00253     fFile->cd();
00254     for(int c=0;c<65;c++) {
00255       fHisto[c] = new TH1F(Form("h%02d",c),
00256                            Form("Pixel %02d",c),
00257                            20500,-500,20000 );
00258     }
00259     if(save) save->cd();
00260   }
00261  
00262 
00263   int ndigits_tot = 0;
00264   int ndigits = 0;
00265 
00266   // Get the RawDigitDataBlock.
00267   TObject* tobj;
00268   TIter    fragiter = mom->FragmentIter();
00269   while( ( tobj = fragiter.Next() ) ) {
00270     RawRecord* rawRec = dynamic_cast<RawRecord*>(tobj);
00271     if(rawRec) {
00272       TIter itr = rawRec->GetRawBlockIter();
00273       RawDataBlock* rawBlk;
00274       while ( ( rawBlk = dynamic_cast<RawDataBlock*>(itr.Next()) ) ) {
00275         RawDigitDataBlock* rddb = dynamic_cast<RawDigitDataBlock*>(rawBlk);
00276         if(rddb) {
00277           
00278           cout << "Ana context: " << rawRec->GetRawHeader()->GetVldContext().AsString() << endl;
00279           PlexHandle plex(rawRec->GetRawHeader()->GetVldContext());
00280 
00281           int npixels = 16;
00282           if(fPsid.GetElecType()==ElecType::kQIE) npixels = 64;
00283           
00284           PlexPixelSpotId pixels[65];
00285           RawChannelId    channels[65];
00286           Double_t        values[65];
00287           for(int i=0;i<npixels;i++){
00288             pixels[i] = fPsid;
00289             pixels[i].SetPixel(i);
00290             channels[i] = plex.GetRawChannelId(pixels[i]);
00291             values[i]=0;
00292           }
00293 
00294           TIter iter = rddb->GetDatumIter();
00295           while( (tobj=iter.Next() ) ) {
00296             
00297             // Process a digit.     
00298             RawDigit*        rawDigit = dynamic_cast<RawDigit *>(tobj);
00299             
00300             if(rawDigit) {
00301               // Try to add it.
00302               // We do it this way to 
00303               // make sure we get all buckets.
00304               ndigits_tot++;
00305               for(int i=0;i<npixels;i++) {
00306                 if(rawDigit->GetChannel().IsSameChannel(channels[i])) {
00307                   ndigits++;
00308                   values[i] += rawDigit->GetADC();
00309                   if(rawDigit->GetChannel().GetElecType()==ElecType::kQIE)
00310                     values[i] -= 50; // Remove digital offset for each bucket.
00311                 }
00312               }    
00313             }
00314 
00315           }
00316 
00317 
00318           MSG("DetSim",Msg::kDebug) << "Got " << ndigits << " digits; total in event = " << ndigits_tot << endl;
00319           // Fill histograms.
00320           for(int i=0;i<npixels;i++) {
00321             fHisto[i]->Fill(values[i]);
00322           }
00323           
00324 
00325 
00326         }       
00327       }
00328     }
00329   } // It's so ugly!
00330   
00331 
00332   return JobCResult::kPassed; // kNoDecision, kFailed, etc.
00333 }