00001 #include "MuonHitSim.h"
00002 #include "Event/SimHeader.h"
00003 #include "Event/SimPmtHit.h"
00004 #include "DataSvc/ICableSvc.h"
00005 #include "Context/ServiceMode.h"
00006 #include "Context/Context.h"
00007 #include "Context/ContextRange.h"
00008
00009 #include "TH1F.h"
00010 #include "TF1.h"
00011 #include "TParameter.h"
00012 #include "TFile.h"
00013 #include "TRandom.h"
00014 #include <string>
00015 #include <vector>
00016
00017 #include "GaudiKernel/IRndmGen.h"
00018 #include "GaudiKernel/IRndmEngine.h"
00019 #include "GaudiKernel/IRndmGenSvc.h"
00020 #include "GaudiKernel/SystemOfUnits.h"
00021 #include "GaudiKernel/PhysicalConstants.h"
00022
00023 using namespace std;
00024 using namespace DayaBay;
00025
00026 MuonHitSim::MuonHitSim(const string& name, ISvcLocator* svcloc)
00027 : GaudiAlgorithm(name, svcloc)
00028 {
00029 declareProperty("MeVtoPE", m_MeVtoPE, "the PE/MeV conversion factor");
00030 declareProperty("CableSvcName", m_cableSvcName = "StaticCableSvc",
00031 "Name of service to map between det, hardware, and electronic IDs");
00032 }
00033
00034 MuonHitSim::~MuonHitSim()
00035 {
00036 }
00037
00038 StatusCode MuonHitSim::initialize()
00039 {
00040 StatusCode sc = this->GaudiAlgorithm::initialize();
00041 if(sc != StatusCode::SUCCESS) return sc;
00042 if(sc.isFailure())
00043 {
00044 error() << "parent class failed to initialize"<< endreq;
00045 return sc;
00046 }
00047
00048
00049 m_cableSvc = svc<ICableSvc>(m_cableSvcName, true);
00050 if(!m_cableSvc)
00051 {
00052 error()<< "Failed to access cable svc: "<< m_cableSvcName<< endreq;
00053 return StatusCode::FAILURE;
00054 }
00055
00056
00057 sc = service("RndmGenSvc", m_ranSvc, true);
00058 if(sc.isFailure())
00059 {
00060 error() << "Unable get RndmGenSvc service"<< m_ranSvc<< endreq;
00061 return sc;
00062 }
00063
00064 string wsFile = getenv("MUONHITSIMROOT");
00065 wsFile.append("/data/Anafile.root");
00066 m_wsFile = new TFile(wsFile.c_str(), "READ");
00067 m_tracPEIWS = (TH2F*)gDirectory->Get("TracPEIWS");
00068 m_tracPEOWS = (TH2F*)gDirectory->Get("TracPEOWS");
00069 m_tracNHitIWS = (TH2F*)gDirectory->Get("TracNHitIWS");
00070 m_tracNHitOWS = (TH2F*)gDirectory->Get("TracNHitOWS");
00071
00072 return StatusCode::SUCCESS;
00073 }
00074
00075 StatusCode MuonHitSim::execute()
00076 {
00077 debug() << "execute() ______________________________ start" << endreq;
00078
00079
00080
00081 SimHeader* simHeader = get<SimHeader>(DayaBay::SimHeaderLocation::Default);
00082
00083
00084 const DayaBay::SimHitHeader* simHitHeader = simHeader->hits();
00085 const DayaBay::SimHitHeader::hc_map& hcMap = simHitHeader->hitCollection();
00086
00087 debug() << "Processing hit collections" <<endreq;
00088
00089 debug() << "SimHeader TimeStamp: " << simHeader->timeStamp() << ", SimHeader Earliest: " << simHeader->earliest() << " , Latest: " <<simHeader->latest() << endreq;
00090
00091
00092
00093 TimeStamp ts = simHeader->timeStamp();
00094
00095 bool existHit = false;
00096 if(!hcMap.empty())
00097 {
00098 existHit = true;
00099 }
00100
00101 m_earliestHit = 300;
00102 m_latestHit = 0;
00103
00104
00105 double depEinLSAD1 = 0;
00106 double depEinGdLSAD1 = 0;
00107 double depEinLSAD2 = 0;
00108 double depEinGdLSAD2 = 0;
00109 double trackinIWS = 0;
00110 double trackinOWS = 0;
00111
00112 const SimUnobservableStatisticsHeader* hstats=simHeader->unobservableStatistics();
00113 map<std::string, SimStatistic>::const_iterator hstatsIter;
00114 for(hstatsIter=hstats->stats().begin(); hstatsIter!=hstats->stats().end(); hstatsIter++)
00115 {
00116 if(hstatsIter->first=="QEDepInLS1")
00117 {
00118 depEinLSAD1=hstatsIter->second.sum();
00119 if(hstatsIter->second.sum()>0)
00120 debug()<<"AD1 qLS1 ="<<hstatsIter->second.sum()<<endreq;
00121 }
00122 if(hstatsIter->first=="QEDepInGdLS1")
00123 {
00124 depEinGdLSAD1=hstatsIter->second.sum();
00125 if(hstatsIter->second.sum()>0)
00126 debug()<<"AD1 qGdLS1 ="<<hstatsIter->second.sum()<<endreq;
00127 }
00128 if(hstatsIter->first=="QEDepInLS2")
00129 {
00130 depEinLSAD2 = hstatsIter->second.sum();
00131 if(hstatsIter->second.sum()>0)
00132 debug()<<"AD2 qLS2 ="<<hstatsIter->second.sum()<<endreq;
00133 }
00134 if(hstatsIter->first=="QEDepInGdLS2")
00135 {
00136 depEinGdLSAD2 = hstatsIter->second.sum();
00137 if(hstatsIter->second.sum()>0)
00138 debug()<<"AD2 qGdLS2 ="<<hstatsIter->second.sum()<<endreq;
00139 }
00140 if(hstatsIter->first=="MuTrackInIWS")
00141 {
00142 trackinIWS = hstatsIter->second.sum();
00143 if(hstatsIter->second.sum()>0)
00144 debug()<<"IWS Track ="<<hstatsIter->second.sum()<<endreq;
00145 }
00146 if(hstatsIter->first=="MuTrackInOWS")
00147 {
00148 trackinOWS = hstatsIter->second.sum();
00149 if(hstatsIter->second.sum()>0)
00150 debug()<<"OWS Track ="<<hstatsIter->second.sum()<<endreq;
00151 }
00152 }
00153
00154 double depEAD1 = 0;
00155 depEAD1 = depEinGdLSAD1 + depEinLSAD1;
00156 debug() << "depE in AD1: " << depEAD1 << endreq;
00157 double depEAD2 = 0;
00158 depEAD2 = depEinGdLSAD2 + depEinLSAD2;
00159 debug() << "depE in AD2: " << depEAD2 << endreq;
00160
00161
00162 debug() << "Track in Inner Water Pool: " << trackinIWS << endreq;
00163 debug() << "Track in Outter Water Pool: " << trackinOWS << endreq;
00164
00165
00166
00167 if(depEAD1+depEAD2+trackinIWS+trackinOWS>0)
00168 {
00169 bool ad1 = false;
00170 bool ad2 = false;
00171 bool iws = false;
00172 bool ows = false;
00173
00175 DayaBay::SimHitHeader::hc_map::const_iterator hcIter;
00176 for(hcIter=hcMap.begin(); hcIter != hcMap.end(); ++hcIter) {
00177
00178 DayaBay::Detector det(hcIter->first);
00179 debug() << "Checking " << det.detName() << " for hits." << endreq;
00180 DayaBay::SimHitCollection* hits = hcIter->second;
00181 if(!hits) return StatusCode::FAILURE;
00182
00183 debug() << "Got hit collection from " << det.detName()<<" (id= "
00184 << det.siteDetPackedData() << ") " << " with "
00185 << hits->collection().size() << " hits." << endreq;
00186
00187 const DayaBay::SimHitCollection::hit_container& htcon = hits->collection();
00188 DayaBay::SimHitCollection::hit_container::const_iterator hcIter,
00189 hcDone = htcon.end();
00190 for (hcIter=htcon.begin(); hcIter != hcDone; ++hcIter) {
00191
00192 DayaBay::DetectorSensor sensDetId((*hcIter)->sensDetId());
00193 debug() << "Requesting pmtData for " << sensDetId << "," << (*hcIter)->hitTime() << endreq;
00194 }
00195
00196 if( det.detName() == "DayaBayAD1") {
00197 ad1 = true;
00198 addPE(det, hits, depEAD1, ts);
00199 } else if( det.detName() == "DayaBayAD2") {
00200 ad2 = true;
00201 addPE(det, hits, depEAD2, ts);
00202 } else if( det.detName() == "DayaBayIWS"){
00203 iws = true;
00204 addWSPE(det, hits, trackinIWS, ts);
00205 } else if( det.detName() == "DayaBayOWS"){
00206 ows = true;
00207 addWSPE(det, hits, trackinOWS, ts);
00208 }
00209 }
00210
00211 if(!ad1 && depEAD1>0)
00212 {
00213 SimHitHeader *shh = const_cast<SimHitHeader*>(simHitHeader);
00214 SimHitCollection *hits = new SimHitCollection;
00215 hits->setHeader(shh);
00216
00217 Detector det("DayaBayAD1");
00218 hits->setDetector(det);
00219
00220 SimHitCollection::hit_container *hc = new SimHitCollection::hit_container;
00221 hits->setCollection(*hc);
00222 addPE(det, hits, depEAD1, ts);
00223 if(hits->collection().empty())
00224 {
00225 delete hits;
00226 }
00227 else
00228 {
00229 shh->addHitCollection(hits);
00230 }
00231 }
00232 if(!ad2 && depEAD2>0)
00233 {
00234 SimHitHeader *shh = const_cast<SimHitHeader*>(simHitHeader);
00235 SimHitCollection *hits = new SimHitCollection;
00236 hits->setHeader(shh);
00237
00238 Detector det("DayaBayAD2");
00239 hits->setDetector(det);
00240
00241 SimHitCollection::hit_container *hc = new SimHitCollection::hit_container;
00242 hits->setCollection(*hc);
00243 addPE(det, hits, depEAD2, ts);
00244 if(hits->collection().empty())
00245 {
00246 delete hits;
00247 }
00248 else
00249 {
00250 shh->addHitCollection(hits);
00251 }
00252 }
00253 if(!iws && trackinIWS>0)
00254 {
00255 SimHitHeader *shh = const_cast<SimHitHeader*>(simHitHeader);
00256 SimHitCollection *hits = new SimHitCollection;
00257 hits->setHeader(shh);
00258
00259 Detector det("DayaBayIWS");
00260 hits->setDetector(det);
00261
00262 SimHitCollection::hit_container *hc = new SimHitCollection::hit_container;
00263 hits->setCollection(*hc);
00264 addWSPE(det, hits, trackinIWS, ts);
00265 if(hits->collection().empty())
00266 {
00267 delete hits;
00268 }
00269 else
00270 {
00271 shh->addHitCollection(hits);
00272 }
00273 }
00274 if(!ows && trackinOWS>0)
00275 {
00276 SimHitHeader *shh = const_cast<SimHitHeader*>(simHitHeader);
00277 SimHitCollection *hits = new SimHitCollection;
00278 hits->setHeader(shh);
00279
00280 Detector det("DayaBayOWS");
00281 hits->setDetector(det);
00282
00283 SimHitCollection::hit_container *hc = new SimHitCollection::hit_container;
00284 hits->setCollection(*hc);
00285 addWSPE(det, hits, trackinOWS, ts);
00286 if(hits->collection().empty())
00287 {
00288 delete hits;
00289 }
00290 else
00291 {
00292 shh->addHitCollection(hits);
00293 }
00294 }
00295
00296
00297
00298
00299
00300
00301
00302
00303
00304
00305
00306
00307
00308
00309
00310
00311
00312
00313
00314
00315
00316
00317
00318
00319
00320
00321
00322
00323
00324
00325
00326
00327
00328
00329
00330
00331
00332
00333 debug() << "Earliest Hit: " << m_earliestHit << ", Latest Hit: " << m_latestHit<< endreq;
00334 if(m_earliestHit <= m_latestHit)
00335 {
00336 TimeStamp tempEarliestStamp = simHeader->timeStamp();
00337 TimeStamp tempearliest(0, m_earliestHit);
00338 tempEarliestStamp.Add(tempearliest);
00339 TimeStamp tempLatestStamp = simHeader->timeStamp();
00340 TimeStamp templatest(0, m_latestHit);
00341 tempLatestStamp.Add(templatest);
00342 if(existHit)
00343 {
00344 if(simHeader->earliest()>tempEarliestStamp)
00345 {
00346 simHeader->setEarliest(tempEarliestStamp);
00347 }
00348 if(simHeader->latest()<tempLatestStamp)
00349 {
00350 simHeader->setLatest(tempLatestStamp);
00351 }
00352 }
00353 else
00354 {
00355 simHeader->setEarliest(tempEarliestStamp);
00356 simHeader->setLatest(tempLatestStamp);
00357 }
00358 }
00359 else if(m_earliestHit != 300 || m_latestHit != 0)
00360 {
00361 error() << "Fail to get earliest hit time or latest hit time!" << endreq;
00362 }
00363
00364
00365
00366
00367
00368
00369
00370 debug() << "After AddPE, Earliest: " << simHeader->earliest() << " , Latest: " <<simHeader->latest() << endreq;
00371
00372
00373
00374
00375 debug() << "Checking SimHitCollection after muon 300ns simulation." << endreq;
00376 for(hcIter=hcMap.begin(); hcIter != hcMap.end(); ++hcIter) {
00377
00378 DayaBay::Detector det(hcIter->first);
00379 debug() << "Checking " << det.detName() << " for hits." << endreq;
00380 DayaBay::SimHitCollection* hits = hcIter->second;
00381 if(!hits) return StatusCode::FAILURE;
00382
00383 debug() << "Got hit collection from " << det.detName()<<" (id= "
00384 << det.siteDetPackedData() << ") " << " with "
00385 << hits->collection().size() << " hits." << endreq;
00386 }
00387 }
00388
00389
00390 debug() << "execute() ______________________________ end" << endreq;
00391
00392 return StatusCode::SUCCESS;
00393 }
00394
00395 StatusCode MuonHitSim::finalize()
00396 {
00397 debug() << "finalize()" << endreq;
00398
00399 if(m_wsFile) {
00400 m_wsFile->Close();
00401 delete m_wsFile;
00402 }
00403
00404 return StatusCode::SUCCESS;
00405 }
00406
00407 void MuonHitSim::addPE(const DayaBay::Detector &det, SimHitCollection* hits,
00408 double depE, const TimeStamp& ts)
00409 {
00410
00411 double totalPE = depE * m_MeVtoPE;
00412 if(totalPE > 198*5000 ) totalPE = 198*5000;
00413 double weightPE = 1;
00414
00415
00416 double PEsum = 0;
00417
00418 Context context;
00419 context.SetSite(det.site());
00420 context.SetDetId(det.detectorId());
00421 context.SetSimFlag(SimFlag::kMC);
00422 context.SetTimeStamp(ts);
00423
00424 ServiceMode svcMode(context, 0);
00425
00426 Rndm::Numbers hittime;
00427 if(hittime.initialize(m_ranSvc, Rndm::Gauss(70, 15)).isSuccess())
00428 {
00429
00430 while(PEsum <= totalPE)
00431 {
00432 const vector<DayaBay::AdPmtSensor>& channelList = m_cableSvc->adPmtSensors(svcMode);
00433 vector<DayaBay::AdPmtSensor>::const_iterator chanIter, chanEnd = channelList.end();
00434 for(chanIter = channelList.begin(); chanIter != chanEnd; chanIter++)
00435 {
00436 DayaBay::AdPmtSensor pmtId = *chanIter;
00437 if(pmtId.detectorId()==det.detectorId())
00438 {
00439 int ring = pmtId.ring();
00440 int column = pmtId.column();
00441 SimPmtHit *hit = new SimPmtHit;
00442 int HTime = (int)hittime();
00443 while(HTime<=0 || HTime>=300)
00444 {
00445 HTime = (int)hittime();
00446 }
00447 if(HTime<m_earliestHit) m_earliestHit = HTime;
00448 if(HTime>m_latestHit) m_latestHit = HTime;
00449 hit->setHc(hits);
00450 hit->setHitTime(HTime);
00451 hit->setWeight(weightPE);
00452 hit->setSensDetId(AdPmtSensor(ring,column,det.site(),det.detectorId()).fullPackedData());
00453 hits->collection().push_back(hit);
00454 PEsum += weightPE;
00455 if(PEsum >= totalPE) break;
00456 }
00457 }
00458 }
00459
00460
00461
00462
00463
00464
00465
00466
00467
00468
00469
00470
00471
00472
00473
00474
00475
00476
00477
00478
00479
00480
00481
00482
00483
00484
00485
00486
00487
00488
00489 }
00490 else
00491 {
00492 error() << "Rndm number for hittime failed to initialize" << endreq;
00493 }
00494 debug() << "Total PMT PE: " << PEsum << endreq;
00495 }
00496
00497
00498 void MuonHitSim::addWSPE(const DayaBay::Detector &det, SimHitCollection* hits,
00499 double track, const TimeStamp& ts)
00500 {
00501 if(track<100) return;
00502 double weightPE = 1;
00503 double PEsum = 0;
00504 Context context;
00505 context.SetSite(det.site());
00506 context.SetDetId(det.detectorId());
00507 context.SetSimFlag(SimFlag::kMC);
00508 context.SetTimeStamp(ts);
00509 ServiceMode svcMode(context, 0);
00510 Rndm::Numbers hittime;
00511 if(hittime.initialize(m_ranSvc, Rndm::Gauss(70, 15)).isSuccess())
00512 {
00513 if( det.detName() == "DayaBayIWS")
00514 {
00515 TH1D* htrackNHitIWS = m_tracNHitIWS->ProjectionY("NHitIWS",(int)(track/100),(int)(track/100));
00516 TH1D* htrackPEIWS = m_tracPEIWS->ProjectionY("yPEIWS",(int)(track/100),(int)(track/100));
00517 int nPMT = 0;
00518 if(htrackNHitIWS->GetEntries()>0)
00519 {
00520 nPMT = (int)htrackNHitIWS->GetRandom();
00521 }
00522 else nPMT = 0;
00523 vector<DayaBay::PoolPmtSensor> channelList = m_cableSvc->poolPmtSensors(svcMode);
00524 int NkickPMT = channelList.size()-nPMT;
00525 debug() << "NHit in IWS:" << nPMT << endreq;
00526 for(int i=1;i<=NkickPMT;i++)
00527 {
00528 Rndm::Numbers kick;
00529 if(kick.initialize(m_ranSvc, Rndm::Flat(0, channelList.size()-1)).isSuccess())
00530 {
00531 channelList.erase(channelList.begin()+(int)kick());
00532 }
00533 else
00534 {
00535 error() << "Rndm number for kick PMT failed to initialize" << endreq;
00536 }
00537 }
00538 vector<DayaBay::PoolPmtSensor>::const_iterator chanIter, chanEnd = channelList.end();
00539 for(chanIter = channelList.begin(); chanIter != chanEnd; chanIter++)
00540 {
00541 DayaBay::PoolPmtSensor pmtId = *chanIter;
00542 if(pmtId.detectorId()==det.detectorId())
00543 {
00544 int wallNumber = pmtId.wallNumber();
00545 int wallSpot = pmtId.wallSpot();
00546 bool inwardFacing = pmtId.inwardFacing();
00547 int nPE = 0;
00548 if(htrackPEIWS->GetEntries()>0)
00549 {
00550 nPE = (int)htrackPEIWS->GetRandom();
00551 }
00552 else nPE = 0;
00553 for(int j=1; j<=nPE; j++)
00554 {
00555 SimPmtHit *hit = new SimPmtHit;
00556 int HTime = (int)hittime();
00557 while(HTime<=0 || HTime>=300)
00558 {
00559 HTime = (int)hittime();
00560 }
00561 if(HTime<m_earliestHit) m_earliestHit = HTime;
00562 if(HTime>m_latestHit) m_latestHit = HTime;
00563 hit->setHc(hits);
00564 hit->setHitTime(HTime);
00565 hit->setWeight(weightPE);
00566 hit->setSensDetId(PoolPmtSensor(wallNumber,wallSpot,inwardFacing,det.site(),det.detectorId()).fullPackedData());
00567 hits->collection().push_back(hit);
00568 PEsum ++;
00569 }
00570 }
00571 }
00572 debug() << "IWS PMT PE: " << PEsum << endreq;
00573 }
00574 if( det.detName() == "DayaBayOWS")
00575 {
00576 TH1D* htrackNHitOWS = m_tracNHitOWS->ProjectionY("NHitOWS",(int)(track/100),(int)(track/100));
00577 TH1D* htrackPEOWS = m_tracPEOWS->ProjectionY("yPEOWS",(int)(track/100),(int)(track/100));
00578 int nPMT = 0;
00579 if(htrackNHitOWS->GetEntries()>0)
00580 {
00581 nPMT = (int)htrackNHitOWS->GetRandom();
00582 }
00583 vector<DayaBay::PoolPmtSensor> channelList = m_cableSvc->poolPmtSensors(svcMode);
00584 int NkickPMT = channelList.size()-nPMT;
00585 debug() << "NHit in OWS:" << nPMT << endreq;
00586 for(int i=1;i<=NkickPMT;i++)
00587 {
00588 Rndm::Numbers kick;
00589 if(kick.initialize(m_ranSvc, Rndm::Flat(0, channelList.size()-1)).isSuccess())
00590 {
00591 channelList.erase(channelList.begin()+(int)kick());
00592 }
00593 else
00594 {
00595 error() << "Rndm number for kick PMT failed to initialize" << endreq;
00596 }
00597 }
00598 vector<DayaBay::PoolPmtSensor>::const_iterator chanIter, chanEnd = channelList.end();
00599 for(chanIter = channelList.begin(); chanIter != chanEnd; chanIter++)
00600 {
00601 DayaBay::PoolPmtSensor pmtId = *chanIter;
00602 if(pmtId.detectorId()==det.detectorId())
00603 {
00604 int wallNumber = pmtId.wallNumber();
00605 int wallSpot = pmtId.wallSpot();
00606 bool inwardFacing = pmtId.inwardFacing();
00607 int nPE = 0;
00608 if(htrackPEOWS->GetEntries()>0)
00609 {
00610 nPE = (int)htrackPEOWS->GetRandom();
00611 }
00612 for(int j=1; j<=nPE; j++)
00613 {
00614 SimPmtHit *hit = new SimPmtHit;
00615 int HTime = (int)hittime();
00616 while(HTime<=0 || HTime>=300)
00617 {
00618 HTime = (int)hittime();
00619 }
00620 if(HTime<m_earliestHit) m_earliestHit = HTime;
00621 if(HTime>m_latestHit) m_latestHit = HTime;
00622 hit->setHc(hits);
00623 hit->setHitTime(HTime);
00624 hit->setWeight(weightPE);
00625 hit->setSensDetId(PoolPmtSensor(wallNumber,wallSpot,inwardFacing,det.site(),det.detectorId()).fullPackedData());
00626 hits->collection().push_back(hit);
00627 PEsum ++;
00628 }
00629 }
00630 }
00631 debug() << "OWS PMT PE: " << PEsum << endreq;
00632 }
00633 }
00634 else
00635 {
00636 error() << "Rndm number for hittime failed to initialize" << endreq;
00637 }
00638 }
00639
