3 #include "StEmcSimulatorMaker.h"
8 #include "StMcEvent/StMcCalorimeterHit.hh"
9 #include "StMcEvent/StMcEmcModuleHitCollection.hh"
10 #include "StMcEvent/StMcEmcHitCollection.hh"
11 #include "StMcEvent/StMcEvent.hh"
12 #include "StMcEvent/StMcTrack.hh"
14 #include "StEvent/StEmcRawHit.h"
15 #include "StEvent/StEmcDetector.h"
16 #include "StEvent/StEmcCollection.h"
17 #include "StEvent/StEvent.h"
19 #include "StMuDSTMaker/COMMON/StMuDst.h"
21 #include "StEmcUtil/database/StBemcTables.h"
22 #include "StEmcUtil/geometry/StEmcGeom.h"
23 #include "StEmcUtil/projection/StEmcPosition.h"
25 #include "StEmcSimpleSimulator.h"
26 #include "StEmcPmtSimulator.h"
32 mEmbeddingMode = (GetMaker(
"emcRaw") || GetMaker(
"Eread"));
33 LOG_INFO <<
"StEmcSimulatorMaker EMBEDDING mode = "<< (Int_t)mEmbeddingMode <<endm;
36 for(
int i=0; i<MAXDETBARREL; i++) {
37 mMakeFullDetector[i] =
false;
38 mCheckStatus[i] =
true;
39 mDoZeroSuppression[i] = (mEmbeddingMode) ?
false:
true;
40 mPedestalCut[i] = 1.5;
41 mCalibOffset[i] = 0.0;
42 mCalibSpread[i] = 0.0;
43 mMaxAdcSpread[i] = 0.0;
53 mDoZeroSuppression[0] =
false;
56 mMakeFullDetector[0] =
true;
59 mEmcCollection = NULL;
61 for(
int i=0; i<MAXDETBARREL; i++) {
63 mGeom[i] = StEmcGeom::instance(i+1);
65 LOG_FATAL <<
"Geometry for detector "<<i+1<<
" undefined" << endm;
74 mSimulatorMode[BTOW-1] = StEmcVirtualSimulator::kPrimarySecondaryFullMode;
75 mSimulatorMode[BPRS-1] = StEmcVirtualSimulator::kPrimarySecondaryFullMode;
76 mSimulatorMode[BSMDE-1] = StEmcVirtualSimulator::kSimpleMode;
77 mSimulatorMode[BSMDP-1] = StEmcVirtualSimulator::kSimpleMode;
79 mIsBFC = GetParentChain()->InheritsFrom(
"StBFChain");
80 LOG_INFO <<
"Is StEmcSimulatorMaker in BFC? " << mIsBFC << endm;
83 StEmcSimulatorMaker::~StEmcSimulatorMaker() {
85 delete mSimulator[BTOW-1];
86 delete mSimulator[BPRS-1];
87 delete mSimulator[BSMDE-1];
88 delete mSimulator[BSMDP-1];
96 mSimulator[BTOW-1] =
new StEmcPmtSimulator(kBarrelEmcTowerId, mSimulatorMode[BTOW-1]);
97 mSimulator[BPRS-1] =
new StEmcPmtSimulator(kBarrelEmcPreShowerId, mSimulatorMode[BPRS-1]);
101 for(
int det=BTOW; det<=BSMDP; det++) {
102 mSimulator[det-1]->setEmbeddingMode(mEmbeddingMode);
103 mSimulator[det-1]->setTables(mTables);
104 mSimulator[det-1]->setCalibScale(1.0 + mCalibOffset[det-1]);
105 mSimulator[det-1]->setCalibSpread(mCalibSpread[det-1]);
106 mSimulator[det-1]->setMaximumAdc(mMaxAdc[det-1]);
107 mSimulator[det-1]->setMaximumAdcSpread(mMaxAdcSpread[det-1]);
118 mMcEvent =
dynamic_cast<StMcEvent*
>(GetDataSet(
"StMcEvent"));
120 LOG_ERROR <<
"couldn't find StMcEvent for this event" << endm;
124 mEmcMcHits[0] = mMcEvent->bemcHitCollection();
125 mEmcMcHits[1] = mMcEvent->bprsHitCollection();
126 mEmcMcHits[2] = mMcEvent->bsmdeHitCollection();
127 mEmcMcHits[3] = mMcEvent->bsmdpHitCollection();
130 vector<StMcTrack*> McTracks = mMcEvent->tracks();
131 for(
unsigned int i = 0; i < McTracks.size(); ++i) {
133 if(track->stopVertex())
continue;
134 makeCrossTalk(track);
138 int maxChannels[4] = {4800, 4800, 18000, 18000};
139 std::vector<int> hasHit;
140 std::vector<int>::const_iterator iter;
142 int softId, module, eta, sub;
143 for(
int det=BTOW; det<=BSMDP; det++) {
144 if(mMakeFullDetector[det-1] == 0)
continue;
147 for(
unsigned int mod=1; mod<=mEmcMcHits[det-1]->numberOfModules(); mod++) {
149 const vector<StMcCalorimeterHit*> hits = module->detectorHits();
150 for(
unsigned long i=0; i<hits.size(); i++) {
151 mGeom[det-1]->getId(hits[i]->module(), hits[i]->eta(), hits[i]->sub(), softId);
152 hasHit.push_back(softId);
157 std::sort(hasHit.begin(), hasHit.end());
158 iter = hasHit.begin();
160 if(hasHit.size()) nextHitId = *iter;
165 for(
int softId=1; softId<=maxChannels[det-1]; softId++) {
166 if( softId != nextHitId ) {
167 mGeom[det-1]->
getBin(softId,module,eta,sub);
168 mcHit->setModule(module);
172 mcHit->setParentTrack(NULL);
175 mEmcMcHits[det-1]->module(module)->detectorHits().push_back(mcHit);
180 if(iter != hasHit.end()) nextHitId = *iter;
189 for(
int det=BTOW; det<=BSMDP; det++) {
190 LOG_DEBUG << *(mEmcMcHits[det-1]) << endm;
191 for(
unsigned int mod=1; mod<=mEmcMcHits[det-1]->numberOfModules(); mod++) {
193 const vector<StMcCalorimeterHit*> hits = module->detectorHits();
194 for(
unsigned long i=0; i<hits.size(); i++) {
195 int softId; mGeom[det-1]->getId(hits[i]->module(), hits[i]->eta(), hits[i]->sub(), softId);
196 LOG_DEBUG <<
"softId:" << softId <<
" mod:" << hits[i]->module() <<
" eta:" << hits[i]->eta() <<
" sub:" << hits[i]->sub() <<
" dE:" << hits[i]->dE() << endm;
208 void StEmcSimulatorMaker::makeRawHits() {
212 if (mEmbeddingMode) {
220 mEmcCollection =
event->emcCollection();
221 if (!mEmcCollection) {
223 event->setEmcCollection(mEmcCollection);
231 for (
int det = BTOW;det <= BSMDP;det++) {
232 StDetectorId detectorId = kUnknownId;
234 case BTOW: detectorId = kBarrelEmcTowerId;
break;
235 case BPRS: detectorId = kBarrelEmcPreShowerId;
break;
236 case BSMDE: detectorId = kBarrelSmdEtaStripId;
break;
237 case BSMDP: detectorId = kBarrelSmdPhiStripId;
break;
240 mEmcCollection->setDetector(detector);
242 for (
unsigned int mod = 1;mod <= mEmcMcHits[det-1]->numberOfModules();mod++) {
244 const vector<StMcCalorimeterHit*> hits = module->detectorHits();
245 for (
unsigned long i = 0;i < hits.size();i++) {
247 mGeom[det-1]->getId(hits[i]->module(), hits[i]->eta(), hits[i]->sub(), softId);
248 LOG_DEBUG <<
"-----------------------------------------------------------------------------------------------" << endm;
250 if (mCheckStatus[det-1]) {
251 if (mTables->
status(det, softId) != 1) {
252 LOG_DEBUG << Form(
"det=%2d softId=%5d -- removing hit b/c DB status!=1", detectorId, softId) << endm;
257 StEmcRawHit *rawHit = mSimulator[det-1]->makeRawHit(hits[i]);
258 rawHit->setCalibrationType(0);
261 if (mDoZeroSuppression[det-1]) {
262 float pedMean = mTables->pedestal(det, softId);
263 float pedRMS = mTables->pedestalRMS(det, softId);
264 if ((rawHit->adc() - pedMean) < (mPedestalCut[det-1] * pedRMS)) {
265 LOG_DEBUG <<
"removing hit b/c it failed pedestal cut" << endm;
272 if (mIsBFC && !mEmbeddingMode) {
273 rawHit->setEnergy(hits[i]->dE());
276 detector->addHit(rawHit);
289 for (
int det = BSMDE;det <= BSMDP;++det) {
291 if(det == BSMDP)
continue;
293 vector<StMcCalorimeterHit*> trackHits;
295 case BSMDE: trackHits = track->bsmdeHits();
break;
296 case BSMDP: trackHits = track->bsmdpHits();
break;
298 LOG_DEBUG <<
" this track has BSMD of size = " << trackHits.size() << endm;
303 double highEnergy = -1;
304 int highElement = -1;
305 for (
unsigned long j = 0;j < trackHits.size();++j) {
306 double energy = trackHits.at(j)->dE();
307 if ((energy > highEnergy) || (highEnergy < 0)) {
312 if (highElement != -1) {
318 const Int_t numCross = 5;
319 int softIds[2 * numCross + 1] = {0};
328 for (Int_t i = 0;i < (2*numCross + 1);i++) {
329 softIds[i] = mPosition->
getNextId(det, highHit->module(), highHit->eta(), highHit->sub(), (det == BSMDE) ? (i - numCross) : 0, (det == BSMDP) ? (i - numCross) : 0);
336 mGeom[det - 1]->getId(highHit->module(), highHit->eta(), highHit->sub(), softIds[numCross]);
337 modHigh = highHit->module();
338 mGeom[det - 1]->getEta(softIds[numCross], etaHigh);
343 double totalenergyb = 0;
344 double totalenergya = 0;
345 for (
unsigned int mod = 1;mod <= mEmcMcHits[det- 1]->numberOfModules();++mod) {
347 const vector<StMcCalorimeterHit*> &detectorHits = module->detectorHits();
348 for (
unsigned long k = 0;k < detectorHits.size();++k) {
350 mGeom[det - 1]->getId(detectorHits[k]->module(), detectorHits[k]->eta(), detectorHits[k]->sub(), softTemp);
352 for (Int_t i = 0;i < (2*numCross + 1);i++) {
353 if (softTemp == softIds[i]) {
354 detectorNextHits[i] = detectorHits[k];
355 totalenergyb += detectorHits[k]->dE();
362 const float crossTalk = mCrossTalk[det - 1] * (1 - fabs(etaHigh) );
364 double leakE1 = 0, leakE2 = 0;
366 for (
int i = 0;i <= numCross;i++) {
367 double energy1 = detectorNextHits[numCross + i] ? detectorNextHits[numCross + i]->dE() : 0;
368 double energy2 = detectorNextHits[numCross - i] ? detectorNextHits[numCross - i]->dE() : 0;
369 {LOG_DEBUG <<
"In: softId1 = " << softIds[numCross + i] <<
", softId2 = " << softIds[numCross - i] << endm;}
370 {LOG_DEBUG <<
"In: energy1 = " << energy1 <<
", energy2 = " << energy2 << endm;}
371 {LOG_DEBUG <<
"In: leakE1 = " << leakE1 <<
", leakE2 = " << leakE2 << endm;}
380 if ((i < numCross) && softIds[numCross + i] && softIds[numCross + (i + 1)]) {
382 const double energyNext1 = detectorNextHits[numCross + (i + 1)] ? detectorNextHits[numCross + (i + 1)]->dE() : 0;
384 leakE1 = crossTalk * (energy1 - energyNext1);
389 if ((i < numCross) && softIds[numCross - i] && softIds[numCross - (i + 1)]) {
391 const double energyNext2 = detectorNextHits[numCross - (i + 1)] ? detectorNextHits[numCross - (i + 1)]->dE() : 0;
393 leakE2 = crossTalk * (energy2 - energyNext2);
404 {LOG_DEBUG <<
"Out: energy1 = " << energy1 <<
", energy2 = " << energy2 << endm;}
405 {LOG_DEBUG <<
"Out: leakE1 = " << leakE1 <<
", leakE2 = " << leakE2 << endm;}
407 int module, eta, sub;
408 if (detectorNextHits[numCross + i]) {
409 detectorNextHits[numCross + i]->setdE(energy1);
411 if (softIds[numCross + i]) {
416 mGeom[det - 1]->
getBin(softIds[numCross + i], module, eta, sub);
417 tempHit->setModule(module);
418 tempHit->setEta(eta);
419 tempHit->setSub(sub);
420 tempHit->setParentTrack(NULL);
421 tempHit->setdE(energy1);
422 detectorNextHits[numCross + i] = tempHit;
423 StMcEmcHitCollection::EAddHit returnCode = mEmcMcHits[det - 1]->addHit(tempHit);
424 if (returnCode == StMcEmcHitCollection::kNew) {
430 if (detectorNextHits[numCross - i]) {
431 detectorNextHits[numCross - i]->setdE(energy2);
433 if (softIds[numCross - i]) {
434 mGeom[det - 1]->
getBin(softIds[numCross - i], module, eta, sub);
435 tempHit->setModule(module);
436 tempHit->setEta(eta);
437 tempHit->setSub(sub);
438 tempHit->setParentTrack(NULL);
439 tempHit->setdE(energy2);
440 detectorNextHits[numCross - i] = tempHit;
441 StMcEmcHitCollection::EAddHit returnCode = mEmcMcHits[det - 1]->addHit(tempHit);
442 if (returnCode == StMcEmcHitCollection::kNew) {
448 totalenergya += detectorNextHits[numCross + i] ? detectorNextHits[numCross + i]->dE() : 0;
449 if (i != 0) totalenergya += detectorNextHits[numCross - i] ? detectorNextHits[numCross - i]->dE() : 0;
451 {LOG_DEBUG <<
"total energy before is " << totalenergyb <<
" and total energy after immediately, is " << totalenergya << endm;}
virtual void AddData(TDataSet *data, const char *dir=".data")
User methods.
int status(int det, int softId, const char *option="") const
Int_t getNextId(Int_t det, Int_t m, Int_t e, Int_t s, Int_t nEta, Int_t nPhi) const
Return neighbor id (works for all detectors 1=bemc, 2=bprs, 3=bsmde, 4=bsmdp)
Monte Carlo Track class All information on a simulated track is stored in this class: kinematics...
void loadTables(StMaker *anyMaker)
load tables.
virtual void Clear(const char *)
resets the pointer to the current StMcEvent
virtual Int_t Init()
creates the detector simulators and sets their properties
Int_t getBin(const Float_t phi, const Float_t eta, Int_t &m, Int_t &e, Int_t &s) const
static void setEmcCollection(StEmcCollection *emc_coll)
set pointer to current StEmcCollection
Event data structure to hold all information from a Monte Carlo simulation. This class is the interfa...