00001
00002
00003
00004
00005
00006
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];
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
00075
00076 static Registry r;
00077
00078
00079 std::string name = this->GetName();
00080 name += ".config.default";
00081 r.SetName(name.c_str());
00082
00083
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
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();
00136
00137
00138
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
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
00154 fContext = VldContext(Detector::CharToEnum((char)fDetector),
00155 SimFlag::kMC,
00156 fPulseTime);
00157
00158
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");
00180
00181
00182 mom -> AdoptFragment(simsnarl);
00183
00184
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
00193 gRandom->SetSeed(fCurrentSnarl);
00194 Int_t npe = gRandom->Poisson(fPulseHeight);
00195
00196
00197 for(int i=0;i<npe;i++) {
00198
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
00220
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 , int snarl, int )
00237 {
00238 fCurrentSnarl = snarl;
00239 return this->Get();
00240 }
00241
00242
00243
00244 JobCResult LightInjector::Ana(const MomNavigator* mom)
00245 {
00246
00247
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
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
00298 RawDigit* rawDigit = dynamic_cast<RawDigit *>(tobj);
00299
00300 if(rawDigit) {
00301
00302
00303
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;
00311 }
00312 }
00313 }
00314
00315 }
00316
00317
00318 MSG("DetSim",Msg::kDebug) << "Got " << ndigits << " digits; total in event = " << ndigits_tot << endl;
00319
00320 for(int i=0;i<npixels;i++) {
00321 fHisto[i]->Fill(values[i]);
00322 }
00323
00324
00325
00326 }
00327 }
00328 }
00329 }
00330
00331
00332 return JobCResult::kPassed;
00333 }