StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
StEEmcSlowMaker.cxx
1 // *-- Author : Hal Spinka
2 //
3 // $Id: StEEmcSlowMaker.cxx,v 2.11 2010/09/07 22:24:52 stevens4 Exp $
4 
5 #include <TFile.h>
6 #include <TH2.h>
7 #include <TRandom.h>
8 #include <StMessMgr.h>
9 
10 #include "StEventTypes.h"
11 #include "StMuDSTMaker/COMMON/StMuTypes.hh"
12 
13 #include "StEEmcUtil/database/EEmcDbItem.h"
14 #include "StEEmcUtil/database/StEEmcDb.h"
15 
16 #include "StEEmcFastMaker.h"
17 #include "StEEmcSlowMaker.h"
18 
19 ClassImp(StEEmcSlowMaker)
20 
21 //________________________________________________
22 StEEmcSlowMaker::StEEmcSlowMaker(const Char_t *name, const Char_t*)
23 : StMaker(name) {
24  mMip2ene = getMipdEdx()*0.7; // This is the SMD thickness of 7 mm
25  // times the minimum ionizing energy loss of
26  // 1.998 MeV/cm from the PDG book
27  mSig1pe = 0.85; // from info from S. Vigdor on MAPMT test results
28  //set different thicknesses for pre and post layers (found in geometry debug by Jason et. al.)
29  mPmip2ene[0] = getMipdEdx()*0.475; // The pre-shower tiles are only 4.75 mm thick.
30  mPmip2ene[1] = getMipdEdx()*0.475; // The pre-shower tiles are only 4.75 mm thick.
31  mPmip2ene[2] = getMipdEdx()*0.5; // The post-shower tiles are only 5 mm thick.
32  mPmip2pe = 2.6*1.5; // 2.6 mip/tower scint * 1.5 light yield
33  // in pre- and post-shower elements
34  // loop to init mMip2pe[] - this will eventually need to be
35  // replaced by a table of measured values
36  for (Int_t i = 0;i < MaxSmdStrips;i++) {
37  // This is currently a trapezoidal-parameterization of
38  // the light curve measurements
39  mMip2pe[i] = avgNumPePerMip(i);
40  }
41 
42  mEeDb=0;
43  mNInpEve=0;
44  memset(mHist,0,sizeof(mHist));
45 
47  mEnableTower=true;
48  mEnableSMD=true;
49  mEnablePrePost=true;
50 
52  mAddPed = true;
54  mSmearPed = true;
56  mDropBad = false;
58  mOverwrite = true;
59 
60  mIsEmbeddingMode = false;
62  mSource = kMuDst;
63 
65  mTruncatePedSmear = 3;
66 
69  setDoLightYield(true);
70 
77  setRelativeLightYield( (4.0/4.75) * 1.68, (4.0/4.75) * 1.68, (4.0/5.0)* 0.94);
78 
80  mSamplingFraction = StEEmcFastMaker::getSamplingFraction();
81  mSamplingFractionUser = mSamplingFraction;
82 
84  for (Int_t i = 0;i < kEEmcNumEtas;i++) {
85  mTowerGains[i] = StEEmcFastMaker::getTowerGains()[i] * mSamplingFraction;
86  }
87  mPrepostGains = StEEmcFastMaker::getPreshowerGain();
88  mSmdGains = StEEmcFastMaker::getSmdGain();
89  mMaxAdc = StEEmcFastMaker::getMaxAdc();
90 
91  if ( (mPrepostGains-23000.0)>1. ||
92  (mSmdGains-23000.0)>1. ||
93  (mTowerGains[0]-18.9654)>0.01 )
94  {
95  LOG_WARN << "--------------------------------------------" << endm;
96  LOG_WARN << "You changed the gains in the fast simulator." << endm;
97  LOG_WARN << "You are on your own." << endm;
98  LOG_WARN << " -- the eemc software team" << endm;
99  LOG_WARN << endm;
100  LOG_WARN << "mSmdGains=" << mSmdGains << endm;
101  LOG_WARN << "mPrepostGains="<<mPrepostGains<<endm;
102  for ( Int_t ii=0;ii<12;ii++ )
103  LOG_WARN << "mTowerGains[etabin="<<ii<<"]="<<mTowerGains[ii]<<endm;
104  //assert(2+2==5);
105  }
106 
107  // initialize tower gain factors to 1
108  for (Int_t sec = 0;sec < kEEmcNumSectors;sec++) {
109  for (Int_t sub = 0;sub < kEEmcNumSubSectors;sub++) {
110  for (Int_t eta = 0;eta < kEEmcNumEtas;eta++) {
111  mTowerGainFact[sec][sub][eta] = 1.0;
112  }
113  }
114  }
115 
116  // initialize SMD gain factors to 1
117  for (Int_t sec = 0;sec < kEEmcNumSectors;sec++) {
118  for (Int_t uv = 0;uv < kEEmcNumSmdUVs;uv++) {
119  for (Int_t strip = 0;strip < kEEmcNumStrips;strip++) {
120  mSmdGainFact[sec][uv][strip] = 1.0;
121  }
122  }
123  }
124 
125  // flag set if used in BFC to change defaults in Init()
126  mIsBFC = GetParentChain()->InheritsFrom("StBFChain");
127 }
128 
129 //________________________________________________
130 StEEmcSlowMaker::~StEEmcSlowMaker() {
131 }
132 
133 //________________________________________________
134 Float_t StEEmcSlowMaker::avgNumPePerMip(Int_t stripID) {
135 //
136 // A parameterization of the average number of photoelectrons
137 // per mip for a given SMD strip. See elog 457.
138 //
139  Float_t ya=0,yb=0,xa=1,xb=2;
140  if ( stripID<1) {
141  ;
142  } else if (stripID<20) {
143  xa=1; ya=2.;
144  xb=20; yb=4.;
145  } else if (stripID<250) {
146  xa=20; ya=4.;
147  xb=250; yb=6.;
148  } else {
149  xa=250; ya=6.;
150  xb=290; yb=9.;
151  }
152  const Float_t y = ya + (yb - ya) / (xb - xa) * (stripID - xa);
153  return y;
154 }
155 
156 //________________________________________________
157 Int_t StEEmcSlowMaker::Init() {
158  LOG_INFO << "mIsEmbeddingMode=" << mIsEmbeddingMode << ", mIsBFC=" << mIsBFC << endm;
159 
160  if (mIsEmbeddingMode) {
161  setDropBad(1); // force bad channels to be dropped in db
162  setAddPed(0); // disable pedestal addition
163  setSmearPed(0); // disable pedestal smearing
164  setOverwrite(1); // overwrites StEvent values
165  setSource("StEvent");
166  } else if (mIsBFC) { //set defaults for BFC with no embedding
167  setAddPed(1); // add pedestals
168  setSmearPed(1); // smear pedestals
169  setDropBad(1); // force bad channels to be dropped in db
170  setOverwrite(1); // overwrite ADC values
171  setSource("StEvent");
172  } else { //set defaults for running in analysis chain
173  if (mSource == kMuDst) disableTower(); //don't change tower in analysis chain since only ADC is stored in MuDst
174  }
175 
176  //print out full configuration in log file
177  LOG_INFO << "mAddPed=" << mAddPed << ", mSmearPed=" << mSmearPed << ", mDropBad=" << mDropBad << ", mOverwrite=" << mOverwrite
178  << ", mEnableTower=" << mEnableTower << ", mEnableSMD=" << mEnableSMD << ", mEnablePrePost=" << mEnablePrePost << endm;
179 
180  mEeDb = (StEEmcDb*)GetDataSet("StEEmcDb");
181  if (!mEeDb) {
182  LOG_ERROR << "Cannot find EEMC database StEEmcDb" << endm;
183  }
184  InitHisto();
185  if ( mSmearPed && !mAddPed) {
186  LOG_WARN << "detected mSmearPed=true && mAddPed=false\nwill not work due to ETOW hits storage container in muDst not accepting negative ADC values." << endm;
187  //assert(0);
188  }
189  if ( mSmearPed ) {
190  LOG_INFO << "detected mSmearPed, (be sure peds>N*sig are loaded in DB!), muDst can't store negative ADC for towers, the NUadc will be forced to be non-negative" << endm;
191  }
192  return StMaker::Init();
193 }
194 
195 //________________________________________________
196 void StEEmcSlowMaker::InitHisto() {
197  Char_t tt1[100], tt2[500];
198 
199  sprintf(tt1,"mm");
200  sprintf(tt2,"freq vs. sector ID; sector ID");
201  mHist[12]=new TH1F(tt1,tt2,20,-0.5,19.5);
202 
203  sprintf(tt1,"PreADC");
204  sprintf(tt2,"Revised ADC for Pre/Post Scint.");
205  mHist[0]=new TH1F(tt1,tt2,100,0,200.);
206 
207  sprintf(tt1,"Pre2ADC");
208  sprintf(tt2,"Initial ADC for Pre/Post Scint.");
209  mHist[1]=new TH1F(tt1,tt2,100,0,200.);
210 
211  sprintf(tt1,"SMDADC");
212  sprintf(tt2,"Revised ADC for SMD Strips");
213  mHist[2]=new TH1F(tt1,tt2,100,0,200.);
214 
215  sprintf(tt1,"SMD2ADC");
216  sprintf(tt2,"Initial ADC for SMD Strips");
217  mHist[3]=new TH1F(tt1,tt2,100,0,200.);
218 
219  mHist[4]=new TH1F("hADCtow","Initial ADC for towers",512+32,-32.,512.);
220  mHist[5]=new TH1F("hADCtow2","Finial ADC fot towers",512+32,-32.,512.);
221 
222  //..................
223  sprintf(tt1,"adc");
224  sprintf(tt2," Adc vs strip ID; strip ID; ADC ");
225  mHist[20]=new TH2F(tt1,tt2,30,0,300,100,0,200.);
226 
227  sprintf(tt1,"ADCvsstr");
228  sprintf(tt2," Adc vs, strip ID; strip ID; ADC ");
229  mHist[19]=new TH2F(tt1,tt2,60,0,300,100,0,200.);
230 
231  // add histos to the list (if provided)
232  for (Int_t i = 0; i < maxHist; ++i) {
233  if (mHist[i]) this->AddHist(mHist[i]);
234  }
235 }
236 
237 //________________________________________________
238 Int_t StEEmcSlowMaker::Make() {
239  Int_t result = StMaker::Make();
240  mNInpEve++;
241  LOG_DEBUG << "iEve " << mNInpEve << ", mSource = " << mSource << endm;
242 
243  switch (mSource) {
244  case kMuDst:
246  {
247  if (!GetInputDS("MuDst")) {
248  LOG_DEBUG<<"::Make() MuDst missing"<<endm;
249  return kStWarn;
250  }
251 
252  StMuEmcCollection* emc = StMuDst::muEmcCollection();
253  if (!emc) {
254  LOG_DEBUG<<"::Make() StMuEmcCollection missing"<<endm;
255  return kStWarn;
256  }
257 
259  if (mEnableTower && (result == kStOk)) result = MakeTower(emc);
260 
262  if (mEnablePrePost && (result == kStOk)) result = MakePrePost(emc);
263 
265  if (mEnableSMD && (result == kStOk)) result = MakeSMD(emc);
266  }
267  break;
268 
269  case kStEvent:
271  {
272 
273  StEmcCollection *emc =0;
274  if(mIsEmbeddingMode) {
275  StEEmcFastMaker *fast = (StEEmcFastMaker*)GetMakerInheritsFrom("StEEmcFastMaker");
276  if(fast==0) {
277  LOG_WARN << GetName() << "::Make() no EEmcFastSim in the chain, ignore Endcap"<< endm;
278  return kStWarn;
279  }
280  emc = fast->GetLocalEmcCollection();
281  } else { // it is not Embedding mode
282 
283  StEvent* event = (StEvent*)GetInputDS("StEvent");
284  if (!event) {
285  LOG_WARN << GetName() << "::Make() no StEvent"<< endm;
286  return kStWarn;
287  }
288  emc = event->emcCollection();
289  }
290  // now emc-collection should be accessible one way or another
291 
292  if (!emc) {
293  LOG_WARN << GetName() << "::Make() no emcCollection()" << endm;
294  return kStWarn;
295  }
296 
298  if (mEnableTower && (result == kStOk)) result = MakeTower(emc);
299 
301  if (mEnablePrePost && (result == kStOk)) result = MakePrePost(emc);
302 
304  if (mEnableSMD && (result == kStOk)) result = MakeSMD(emc);
305  }
306  break;
307 
308  default:
309  LOG_ERROR<< "Unknown source type " << mSource << " for this event" << endm;
310  break;
311  }
312  return result;
313 }
314 
315 //________________________________________________
316 //
317 // Some notes on differences between using StEvent and StMuDst
318 //
319 // Roundoff errors
320 //
321 // Since the MuDst only stores ADC values for the towers, it will be
322 // subject to roundoff errors. This is most noticable for towers with
323 // small ADC response. Keep in mind that (int)2.7 rounds down to 2.
324 //
325 // To partially offset this, we will always assume a starting ADC value
326 // of ADC+0.5 for the MuDst.
327 //
328 
329 Int_t StEEmcSlowMaker::MakeTower(StMuEmcCollection *emc) {
330 
331  /**************************************************************
332  *
333  * Performs pedestal addition and smearing, computes ADC based
334  * on GEANT energy deposit and database gains, and corrects
335  * tower ADC for the relative brightness of the pre- and post-
336  * shower layers.
337  *
338  **************************************************************
339  */
340 
342  Float_t prepost_adc[kEEmcNumSectors][kEEmcNumSubSectors][kEEmcNumEtas];
343  memset( prepost_adc, 0, sizeof(prepost_adc) );
344 
348  if (mDoLightYield) for (Int_t i = 0;i < emc->getNEndcapPrsHits();i++) {
350  Int_t sec,sub,eta,pre;
351  StMuEmcHit *hit=emc->getEndcapPrsHit(i,sec,sub,eta,pre);
352  // range check on returned values
353  if (!( sec >= 1 && sec <= 12) || !(sub >= 1 && sub <= 5) || !(eta >= 1 && eta <= 12) || !(pre >= 1 && pre <= 3)) {
354  LOG_ERROR << "Indexing errors detected: EPRS hit " << i << ", sec = " << sec << ", sub = " << sub << ", eta = " << eta << ", pre = " << pre << endm;
355  setZeroAdc(emc);
356  return kStErr;
357  }
360  const Float_t edeposit = hit ? hit->getEnergy() : 0;
361 
362  // divide by sampling fraction, multiply by _tower_ gain
363  // and brightness factor, and add to the adc sum for
364  // this tower. Since GEANT already accounts for the
365  // difference in thickness between pre/post and normal
366  // layers, the factor of 0.8 is introduced to prevent us
367  // from double correcting.
368  const EEmcDbItem *tower = mEeDb ? mEeDb->getTile(sec,sub-1+'A', eta, 'T') : 0;
369  if (!tower) {
370  LOG_ERROR << "Cannot find DB entry for ETOW: sec = " << sec << ", sub = " << sub << ", eta = " << eta << endm;
371  continue;
372  }
373  Float_t myadc = edeposit / mSamplingFractionUser * tower->gain;
374  myadc *= ( mRelativeLightYield[pre-1] - 1.0 );
375 
376  prepost_adc[sec-1][sub-1][eta-1] += myadc;
377  }
378 
382  for (Int_t i = 0;i < emc->getNEndcapTowerADC();i++) {
383  Int_t sec,sub,eta,adc;
385  emc->getEndcapTowerADC(i,adc,sec,sub,eta);
386  // range check on returned values
387  if (!(sec >= 1 && sec <= 12) || !(sub >= 1 && sub <= 5) || !(eta >= 1 && eta <= 12)) {
388  LOG_ERROR << "Indexing errors detected: ETOW hit " << i << ", sec = " << sec << ", sub = " << sub << ", eta = " << eta << endm;
389  setZeroAdc(emc);
390  return kStErr;
391  }
392 
393  const Int_t old = adc;
394  if (mHist[4]) mHist[4]->Fill( old );
395 
397  const EEmcDbItem *tower = mEeDb ? mEeDb->getTile(sec,sub-1+'A', eta, 'T') : 0;
398  if (!tower) {
399  LOG_ERROR << "Cannot find DB entry for ETOW: sec = " << sec << ", sub = " << sub << ", eta = " << eta << endm;
400  continue;
401  }
402 
403  if (mDropBad && (tower->fail || !checkDBped(tower))) {
404  //zero out values if tower has fail bit set or ped is bad
405  adc = 0;
406  } else {
408  Float_t myadc=(Float_t)adc + 0.5;
411  myadc *= tower->gain / mTowerGains[eta-1] * mSamplingFraction / mSamplingFractionUser;
413  myadc += prepost_adc[sec-1][sub-1][eta-1];
415  myadc *= ( mTowerGainFact[sec-1][sub-1][eta-1] );
416  Float_t ped = tower->ped;
418  if (mAddPed) {
419  if ( mSmearPed ) ped += getPedSmear( tower->sigPed );
420  myadc += ped;
421  }
423  adc = (Int_t)myadc;
424  if (mHist[5]) mHist[5]->Fill(adc);
425  if (mDropBad && (tower->gain <= 0) && mAddPed) {
426  //still add ped to tower if gain is bad since trigger emulator expects peds
427  adc = ped;
428  }
430  if (adc < 0) adc = 0;
431  if (adc > mMaxAdc) adc = mMaxAdc;
432  }
433  if (mOverwrite) {
434  if (old) {
435  LOG_DEBUG << "overwriting tower=" << tower->name << " old adc=" << old << " new adc=" << adc << endm;
436  }
442  emc->setTowerADC( i+1, adc, eemc );
443  }
444  }
445  return kStOk;
446 }
447 
448 //________________________________________________
449 Int_t StEEmcSlowMaker::MakePrePost(StMuEmcCollection *emc) {
450  for (Int_t i = 0;i < emc->getNEndcapPrsHits();i++) {
451  Int_t pre,sec,eta,sub;
453  StMuEmcHit *hit = emc->getEndcapPrsHit(i,sec,sub,eta,pre);
454  if (!hit) continue;
455 
456  // range check on returned values
457  if (!( sec >= 1 && sec <= 12) || !(sub >= 1 && sub <= 5) || !(eta >= 1 && eta <= 12) || !(pre >= 1 && pre <= 3)) {
458  LOG_ERROR << "Indexing errors detected: EPRS hit " << i << ", sec = " << sec << ", sub = " << sub << ", eta = " << eta << ", pre = " << pre << endm;
459  setZeroAdc(emc);
460  return kStErr;
461  }
462 
464  if (mEeDb && (sec < mEeDb->getFirstSector() || sec > mEeDb->getLastSector())) continue;
465 
467  const EEmcDbItem *x = mEeDb ? mEeDb->getTile(sec,sub-1+'A', eta, pre-1+'P') : 0;
468  if (!x) {
469  LOG_ERROR << "Cannot find DB entry for EPRS: sec = " << sec << ", sub = " << sub << ", eta = " << eta << ", pre = " << pre << endm;
470  continue;
471  }
472 
473  Int_t NUadc;
474  // Zero out values if tower is marked as bad
475  if (mDropBad && (x->fail || !checkDBped(x) || (x->gain <= 0))) {
476  NUadc=0;
477  } else {
478  const Float_t adc = hit->getAdc();
479  const Float_t energy = hit->getEnergy();
480 
481  Float_t newadc = energy * x->gain;
482  if (adc > 0) {
483  if (mHist[1]) mHist[1]->Fill(adc);
484 
487  const Float_t mipval = energy / mPmip2ene[pre-1];
488  const Float_t avgpe = mipval * mPmip2pe;
489  const Float_t Npe = gRandom->Poisson(avgpe);
490 
491  if (Npe > 0) {
493  const Float_t sigmape = sqrt(Npe) * mSig1pe;
494  Float_t smearedpe = gRandom->Gaus(Npe,sigmape);
495  if (smearedpe < 0) smearedpe = 0;
497  newadc = smearedpe * mMip2ene*(x->gain) / mPmip2pe;
498  }
499  } // non-zero ADC on input
500 
502  Float_t ped = (mAddPed) ? x->ped : 0;
504  if ( mSmearPed ) ped += getPedSmear( x->sigPed );
505 
507  NUadc = (Int_t)(newadc + 0.5 + ped );
508 
510  if (NUadc < 0) NUadc = 0;
511  if (NUadc > mMaxAdc) NUadc = mMaxAdc;
512  if (adc > 0 && mHist[0]) mHist[0]->Fill(NUadc); //??
513  }
514 
518  if (mOverwrite) hit->setAdc(NUadc);
519  }
520  return kStOk;
521 }
522 
523 //________________________________________________
524 Int_t StEEmcSlowMaker::MakeSMD(StMuEmcCollection *emc) {
525  Int_t iuv = 0;
526  for (Char_t uv = 'U';uv <= 'V';uv++) {
527  iuv++;
528  Int_t sec,strip;
529  for (Int_t i = 0;i < emc->getNEndcapSmdHits(uv);i++) {
530  StMuEmcHit *hit=emc->getEndcapSmdHit(uv,i,sec,strip);
531  // range check on returned values
532  if (!(strip-1 >= 0) || !(strip <= 288) || !(iuv-1 >= 0 && iuv-1 < 2) || !(sec > 0 && sec <= MaxSectors)) {
533  LOG_ERROR << "Bad index for ESMD: sec = " << sec << ", iuv = " << iuv << ", strip = " << strip << endm;
534  setZeroAdc(emc);
535  return kStErr;
536  }
537 
538  // tmp, for fasted analysis use only hits from sectors init in DB
539  if (mEeDb && (sec < mEeDb->getFirstSector() || sec > mEeDb->getLastSector())) continue;
540  const EEmcDbItem *x = mEeDb ? mEeDb->getByStrip(sec,uv,strip) : 0;
541  if (!x) {
542  LOG_ERROR << "Cannot find DB entry for ESMD: sec = " << sec << ", uv = " << uv << ", strip = " << strip << endm;
543  continue;
544  }
545 
546  Int_t NUadc;
547  // Zero out values if tower is marked as bad
548  if (mDropBad && (x->fail || !checkDBped(x) || (x->gain <= 0))) {
549  NUadc=0;
550  } else {
551  const Float_t adc = hit->getAdc();
552  const Float_t energy = hit->getEnergy();
553  Float_t newadc = energy * x->gain;
554 
555  if (adc > 0) {
556  if (mHist[12]) mHist[12]->Fill(x->sec);
557  if (mHist[20]) mHist[20]->Fill(x->strip,adc);
558  if (mHist[3]) mHist[3]->Fill(adc);
559 
560  // Addition of Poisson fluctuations and
561  // Gaussian 1 p.e. resolution
562  const Float_t mipval = energy / mMip2ene;
563  const Float_t avgpe = mipval * mMip2pe[strip];
564  const Float_t Npe = gRandom->Poisson(avgpe);
565  if (Npe > 0) {
567  const Float_t sigmape = sqrt(Npe) * mSig1pe;
568  Float_t smearedpe = gRandom->Gaus(Npe, sigmape);
569  if (smearedpe < 0) smearedpe = 0;
571  newadc = smearedpe * mMip2ene * (x->gain) / mMip2pe[strip];
572  }
573  }
574 
575  // Add in factor from gain smearing
576  newadc *= mSmdGainFact[sec-1][iuv-1][strip-1];
577 
578  // Lookup pedestal in database (possibly zero)
579  Float_t ped = (mAddPed) ? x->ped : 0;
581  if ( mSmearPed ) ped += getPedSmear( x->sigPed );
582 
584  NUadc = (Int_t)(newadc + 0.5 + ped );
585 
587  if (NUadc < 0) NUadc = 0;
588  if (NUadc > mMaxAdc) NUadc = mMaxAdc;
589  if (adc > 0 && mHist[19]) mHist[19]->Fill(x->strip,NUadc);
590  if (adc > 0 && mHist[2]) mHist[2]->Fill(NUadc);
591  }
592 
596  if (mOverwrite) hit->setAdc(NUadc);
597  }// loop over 1 plane
598  } // loop over U,V
599  return kStOk;
600 }
601 
602 //________________________________________________
603 void StEEmcSlowMaker::setSource(const Char_t* name) {
604  if (strcmp(name, "MuDst") == 0) {
605  mSource = kMuDst;
606  } else if (strcmp(name, "StEvent") == 0) {
607  mSource = kStEvent;
608  } else {
609  LOG_WARN<<"::setSource()"<<"Source must be \"MuDst\" or \"StEvent\""<<endm;
610  }
611 }
612 
613 //________________________________________________
614 Int_t StEEmcSlowMaker::MakeTower(StEmcCollection* emc) {
615  StEmcDetector *det = emc->detector(kEndcapEmcTowerId);
616  if (!det) {
617  LOG_DEBUG<<"No kEndcapEmcTowerId"<<endm;
618  return kStOk;
619  }
620  StEmcDetector *prepost = emc->detector(kEndcapEmcPreShowerId);
621  if (!prepost) {
622  LOG_WARN<<"No kEndcapEmcPreShowerId, towers are fast simu only."<<endm;
623  return kStOk;
624  }
625 
627  Float_t prepost_adc[kEEmcNumSectors][kEEmcNumSubSectors][kEEmcNumEtas];
628  memset( prepost_adc, 0, sizeof(prepost_adc) );
629 
630  for (UInt_t sec = 1;sec <= det->numberOfModules();sec++) {
631  StSPtrVecEmcRawHit &prepost_hits = prepost->module(sec)->hits();
632  if (mDoLightYield) for (UInt_t ihit = 0;ihit < prepost_hits.size();ihit++) {
633  StEmcRawHit *hit = prepost_hits[ihit];
634  const UInt_t sub = (hit->sub()-1)%5+1;
635  const UInt_t eta = hit->eta();
636  const UInt_t pre = (hit->sub()-1)/5+1;
637  // range check on returned values
638  if (!(sec >= 1 && sec <= 12) || !(sub >= 1 && sub <= 5) || !(eta >= 1 && eta <= 12) || !(pre >= 1 && pre <= 3)) {
639  LOG_ERROR << "Indexing errors detected for EPRS: sec = " << sec << ", sub = " << sub << ", eta = " << eta << ", pre = " << pre << endm;
640  setZeroAdc(emc);
641  return kStErr;
642  }
643 
644  // energy deposit in GeV
645  const Float_t edeposit = hit->energy();
646  // divide by sampling fraction, multiply by _tower_ gain
647  // and brightness factor, and add to the adc sum for
648  // this tower. Since GEANT already accounts for the
649  // difference in thickness between pre/post and normal
650  // layers, the factor of 0.8 is introduced to prevent us
651  // from double correcting.
652  const EEmcDbItem *tower = mEeDb ? mEeDb->getTile(sec,sub-1+'A', eta, 'T') : 0;
653  if (!tower) {
654  LOG_ERROR << "Cannot find DB entry for ETOW: sec = " << sec << ", sub = " << sub << ", eta = " << eta << endm;
655  continue;
656  }
657  Float_t myadc = edeposit / mSamplingFraction * tower->gain;
658  myadc *= ( mRelativeLightYield[pre-1] - 1.0 );
659  prepost_adc[sec-1][sub-1][eta-1] += myadc;
660  }
661  Int_t nTchange = 0; // counts # of towers with ADC (int) values changed
662  StSPtrVecEmcRawHit &tower_hits = det->module(sec)->hits();
663  for (UInt_t ihit = 0;ihit < tower_hits.size();ihit++) {
664  StEmcRawHit *hit = tower_hits[ihit];
665  const UInt_t sub = hit->sub();
666  const UInt_t eta = hit->eta();
667  // range check on returned values
668  if (!(sec >= 1 && sec <= 12) || !(sub >= 1 && sub <= 5) || !(eta >= 1 && eta <= 12)) {
669  LOG_ERROR << "Indexing errors detected for ETOW: sec = " << sec << ", sub = " << sub << ", eta = " << eta << endm;
670  setZeroAdc(emc);
671  return kStErr;
672  }
673 
674  // get DB entry for this tower
675  const EEmcDbItem *tower = mEeDb ? mEeDb->getTile(sec,sub-1+'A', eta, 'T') : 0;
676  if (!tower) {
677  LOG_ERROR << "Cannot find DB entry for ETOW: sec = " << sec << ", sub = " << sub << ", eta = " << eta << endm;
678  continue;
679  }
680  Int_t adc;
681  if (mDropBad && (tower->fail || !checkDBped(tower))) {
682  //zero out values if tower has fail bit set or ped is bad
683  adc = 0;
684  } else {
685  if (mHist[4]) mHist[4]->Fill( hit->adc() );
686  // energy deposit in GeV
687  const Float_t edeposit = hit->energy();
688  // compute adc using gains stored in the databse
689  Float_t myadc = edeposit / mSamplingFraction * tower->gain;
690  // add in brightness correction for the pre and postshower layers
691  myadc += prepost_adc[sec-1][sub-1][eta-1];
693  myadc *= ( mTowerGainFact[sec-1][sub-1][eta-1] );
695  Float_t ped = tower->ped;
696  if (mAddPed) {
697  if ( mSmearPed ) ped += getPedSmear( tower->sigPed );
698  myadc += ped;
699  }
700  adc = (Int_t)myadc;
701  // Zero out values if tower is marked as bad
702  if (mDropBad && (tower->gain <= 0) && mAddPed) {
703  //still add ped to tower if gain is bad since trigger emulator expects peds
704  adc = ped;
705  }
707  if (adc < 0) adc = 0;
708  if (adc > mMaxAdc) adc = mMaxAdc;
709  }
710  if (mHist[5]) mHist[5]->Fill(adc);
711  if (mOverwrite) {
712  if ((Int_t)hit->adc() != adc) nTchange++;
713  LOG_DEBUG <<"overwriting tower=" << tower->name << " old adc=" << hit->adc() << " new adc=" << adc << endm;
714  hit->setAdc( adc );
715  }
716  }
717  LOG_DEBUG << " Etow changed " << nTchange << " ADC values for sector=" << sec << endm;
718  }// loop over sectors
719  return kStOk;
720 }
721 
722 //________________________________________________
723 Int_t StEEmcSlowMaker::MakePrePost(StEmcCollection* emc) {
724  StEmcDetector* det = emc->detector(kEndcapEmcPreShowerId);
725  if (!det) {
726  LOG_WARN << "No kEndcapEmcPreShowerId" << endm;
727  return kStOk;
728  }
729 
730  for (UInt_t sector =1; sector <= det->numberOfModules(); ++sector) {
731  StSPtrVecEmcRawHit& hits = det->module(sector)->hits();
732  for(UInt_t i = 0; i < hits.size(); ++i) {
733  StEmcRawHit* hit = hits[i];
734  Char_t sub = 'A'+(hit->sub()-1)%5;
735 
736  // Layer: 'P'=preshower1, 'Q'=preshower2, 'R'=postshower
737  Char_t layer = 'P'+(hit->sub()-1)/5;
738  Int_t layerP = (hit->sub()-1)/5;
739 
740  const UInt_t isub = (hit->sub()-1)%5+1;
741  const UInt_t ieta = hit->eta();
742  const UInt_t ipre = (hit->sub()-1)/5+1;
743  // range check on returned values
744  if (!(sector >= 1 && sector <= 12) || !(isub >= 1 && isub <= 5) || !(ieta >= 1 && ieta <= 12) || !(ipre >= 1 && ipre <= 3)) {
745  LOG_ERROR << "Indexing errors detected for EPRS: sec = " << sector << ", sub = " << isub << ", eta = " << ieta << ", pre = " << ipre << endm;
746  setZeroAdc(emc);
747  return kStErr;
748  }
749 
750  // Database ranges: sector=1-12, sub=A-E, eta=1-12, type=T,P-R; Slow method
751  const EEmcDbItem* x = mEeDb ? mEeDb->getTile(sector, sub, hit->eta(), layer) : 0;
752  if (!x) {
753  LOG_ERROR << "Cannot find DB entry for EPRS: sec = " << sector << ", sub = " << sub << ", eta = " << hit->eta() << ", pre = " << layer << endm;
754  continue;
755  }
756 
757  Int_t NUadc;
758  // Zero out values if tower is marked as bad
759  if (mDropBad && (x->fail || !checkDBped(x) || (x->gain <= 0))) {
760  NUadc=0;
761  } else {
762  Float_t adc = hit->adc();
763  Float_t energy = hit->energy();
764  Float_t newadc = energy * x->gain;
765  if(adc > 0) {
766  if (mHist[1]) mHist[1]->Fill(adc);
767 
770  const Float_t mipval = energy / mPmip2ene[layerP];
771  const Float_t avgpe = mipval * mPmip2pe;
772  const Float_t Npe = gRandom->Poisson(avgpe);
773 
774  if (Npe > 0) {
776  const Float_t sigmape = sqrt(Npe) * mSig1pe;
777  Float_t smearedpe = gRandom->Gaus(Npe, sigmape);
778  if (smearedpe < 0) smearedpe = 0;
780  newadc = smearedpe * mMip2ene * x->gain / mPmip2pe;
781  }
782  } // non-zero ADC on input
783 
785  Float_t ped = mAddPed ? x->ped : 0;
786 
788  if ( mSmearPed ) ped += getPedSmear( x->sigPed );
789 
791  NUadc = Int_t(newadc + 0.5 + ped);
792 
794  if ( NUadc < 0 ) NUadc = 0;
795  if ( NUadc > mMaxAdc ) NUadc = mMaxAdc;
796  if (adc > 0 && mHist[0]) mHist[0]->Fill(NUadc); //??
797  }
798 
802  if (mOverwrite) hit->setAdc(NUadc);
803  }
804  }
805  return kStOk;
806 }
807 
808 //________________________________________________
809 Int_t StEEmcSlowMaker::MakeSMD(StEmcCollection* emc) {
810  Int_t iuv = 0;
811  for (Char_t plane = 'U'; plane <= 'V'; ++plane) {
812  iuv++;
813  StEmcDetector* det = 0;
814  switch (plane) {
815  case 'U':
816  det = emc->detector(kEndcapSmdUStripId);
817  if (!det) {
818  LOG_DEBUG << "No kEndcapSmdUStripId" << endm;
819  continue;
820  }
821  break;
822  case 'V':
823  det = emc->detector(kEndcapSmdVStripId);
824  if (!det) {
825  LOG_DEBUG << "No kEndcapSmdVStripId" << endm;
826  continue;
827  }
828  break;
829  }
830 
831  if (!det) continue;
832 
833  for (UInt_t sector = 1; sector <= det->numberOfModules();++sector) {
834  StSPtrVecEmcRawHit& hits = det->module(sector)->hits();
835  for (UInt_t i = 0; i < hits.size(); ++i) {
836  StEmcRawHit* hit = hits[i];
837  Int_t strip = hit->eta();
838 
839  // range check on returned values
840  if (!(strip-1 >= 0) || !(strip <= 288) || !(iuv-1 >= 0 && iuv-1 < 2) || !(sector > 0 && sector <= MaxSectors)) {
841  LOG_ERROR << "Bad index for ESMD: sec = " << sector << ", iuv = " << iuv << ", strip = " << strip << endm;
842  setZeroAdc(emc);
843  return kStErr;
844  }
845 
846  // Database ranges: sector=1-12, plane=U-V, strip=1-288
847  const EEmcDbItem* x = mEeDb ? mEeDb->getByStrip(sector, plane, strip) : 0;
848  if (!x) {
849  LOG_ERROR << "Cannot find DB entry for ESMD: sec = " << sector << ", uv = " << iuv-1+'U' << ", strip = " << strip << endm;
850  continue;
851  }
852 
853  Int_t NUadc;
854  // Zero out values if tower is marked as bad
855  if (mDropBad && (x->fail || !checkDBped(x) || (x->gain <= 0))) {
856  NUadc = 0;
857  } else {
858  const Float_t adc = hit->adc();
859  const Float_t energy = hit->energy();
860  Float_t newadc = energy * x->gain;
861  if (adc > 0) {
862  if (mHist[12]) mHist[12]->Fill(x->sec);
863  if (mHist[2]) mHist[20]->Fill(x->strip, adc);
864  if (mHist[3]) mHist[3]->Fill(adc);
865  // Addition of Poisson fluctuations and
866  // Gaussian 1 p.e. resolution
867  const Float_t mipval = energy / mMip2ene;
868  const Float_t avgpe = mipval * mMip2pe[strip];
869  const Float_t Npe = gRandom->Poisson(avgpe);
870  if (Npe > 0) {
872  const Float_t sigmape = sqrt(Npe) * mSig1pe;
873  Float_t smearedpe = gRandom->Gaus(Npe, sigmape);
874  if (smearedpe < 0) smearedpe = 0;
876  newadc = smearedpe * mMip2ene * x->gain / mMip2pe[strip];
877  }
878  }
879  // Add in factor from gain smearing
880  newadc *= mSmdGainFact[sector-1][iuv-1][strip-1];
882  Float_t ped = mAddPed ? x->ped : 0;
884  if ( mSmearPed ) ped += getPedSmear( x->sigPed );
886  NUadc = Int_t(newadc + 0.5 + ped);
888  if ( NUadc < 0 ) NUadc = 0;
889  if ( NUadc > mMaxAdc ) NUadc = mMaxAdc;
890  if (adc > 0 && mHist[19]) mHist[19]->Fill(x->strip, NUadc);
891  if (adc > 0 && mHist[2]) mHist[2]->Fill(NUadc);
892  }
893 
897  if (mOverwrite) hit->setAdc(NUadc);
898  } // Loop over 1 plane
899  }
900  }
901  return kStOk;
902 }
903 
904 //________________________________________________
905 void StEEmcSlowMaker::setTowerGainSpread(Float_t s, Float_t mean) {
906  LOG_INFO << "setTowerGainSpread(): gain spread: " << s << "; gain mean value: " << mean << endm;
907  // initialize tower gain factors to mean +/- spread
908  for ( Int_t sec=0;sec<kEEmcNumSectors;sec++ ) {
909  for ( Int_t sub=0;sub<kEEmcNumSubSectors;sub++ ) {
910  for ( Int_t eta=0;eta<kEEmcNumEtas;eta++ ) {
911  Float_t f = -1.0E9;
912  while (f <= -mean || f >= 1.0) f = gRandom->Gaus(0, s);
913  mTowerGainFact[sec][sub][eta] = mean + f;
914  }
915  }
916  }
917 }
918 
919 //________________________________________________
920 void StEEmcSlowMaker::setSmdGainSpread(Float_t s) {
921  for (Int_t strip = 0;strip < kEEmcNumStrips;strip++) {
922  setSmdGainSpread(s, strip);
923  }
924 }
925 
926 //________________________________________________
927 void StEEmcSlowMaker::setSmdGainSpread(Float_t s, Int_t strip) {
928  for (Int_t sec = 0;sec < kEEmcNumSectors;sec++) {
929  for (Int_t uv = 0;uv < kEEmcNumSmdUVs;uv++) {
930  setSmdGainSpread(s, sec, uv, strip);
931  }
932  }
933 }
934 
935 //________________________________________________
936 void StEEmcSlowMaker::setSmdGainSpread(Float_t s, Int_t sec, Int_t uv, Int_t strip) {
937  const Float_t mean = 1.0;
938  Float_t f = -1.0E9;
939  while (f <= -mean || f >= 1.0) f = gRandom->Gaus(0, s);
940  mSmdGainFact[sec][uv][strip] = mean + f;
941 }
942 
943 //________________________________________________
944 Bool_t StEEmcSlowMaker::checkDBped(const EEmcDbItem *x) {
945  //check that channels ped and pedSigma are all >= 0
946  return (x && (x->ped >= 0) && (x->sigPed >= 0));
947 }
948 
949 //________________________________________________
950 Float_t StEEmcSlowMaker::getPedSmear(Float_t sigPed) {
951  //get pedestal smearing from gaussian distribution truncated at N*sigma
952  Float_t smear = 0;
953  if (mTruncatePedSmear > 0) {
954  do {
955  smear = gRandom->Gaus(0, 1.0);
956  } while (fabs(smear) > mTruncatePedSmear);
957  }
958  return smear * sigPed;
959 }
960 
961 //________________________________________________
962 void StEEmcSlowMaker::setZeroAdc(StEmcCollection* emc) {
963  StDetectorId detId[4]={kEndcapEmcTowerId,kEndcapEmcPreShowerId,kEndcapSmdUStripId,kEndcapSmdVStripId};
964  StEmcDetector *det=0;
965 
966  for(int i=0; i<4; i++) {
967  det = emc->detector(detId[i]);
968  if (!det) continue;
969  for (UInt_t sec = 1;sec <= det->numberOfModules();sec++) {
970  StSPtrVecEmcRawHit &det_hits = det->module(sec)->hits();
971  for (UInt_t ihit = 0;ihit < det_hits.size();ihit++) {
972  StEmcRawHit *hit = det_hits[ihit];
973  hit->setAdc(0);
974  }
975  }
976  }
977 }
978 
979 //________________________________________________
980 void StEEmcSlowMaker::setZeroAdc(StMuEmcCollection *emc) {
981 
982  for (Int_t i = 0;i < emc->getNEndcapTowerADC();i++)
983  emc->setTowerADC( i+1, 0, eemc );
984 
985  for (Int_t i = 0;i < emc->getNEndcapPrsHits();i++) {
986  Int_t pre,sec,eta,sub;
987  StMuEmcHit *hit = emc->getEndcapPrsHit(i,sec,sub,eta,pre);
988  if (hit) hit->setAdc(0);
989  }
990 
991  for (Char_t uv = 'U';uv <= 'V';uv++) {
992  Int_t sec,strip;
993  for (Int_t i = 0;i < emc->getNEndcapSmdHits(uv);i++) {
994  StMuEmcHit *hit=emc->getEndcapSmdHit(uv,i,sec,strip);
995  if (hit) hit->setAdc(0);
996  }
997  }
998 }
999 
1000 //________________________________________________
1001 Float_t StEEmcSlowMaker::getMipdEdx() {
1002  // Return MIP dE/dx = 1.998 MeV/cm from the PDG book
1003  // used to simulate SMD, Pre, Post ADC response
1004  return 0.001998;
1005 }
1006 
1007 // $Log: StEEmcSlowMaker.cxx,v $
1008 // Revision 2.11 2010/09/07 22:24:52 stevens4
1009 // give access to MIP dE/dx to other makers
1010 //
1011 // Revision 2.10 2010/08/03 02:20:40 stevens4
1012 // final update from peer review
1013 //
1014 // Revision 2.9 2010/07/29 16:12:03 ogrebeny
1015 // Update after the peer review
1016 //
1017 // Revision 2.8 2010/05/03 20:47:22 ogrebeny
1018 // Some code cleanup
1019 //
1020 // Revision 2.7 2010/02/12 23:02:38 ogrebeny
1021 // By the request of the photon group, added an option to shift EEMC gains in the slow simulator.
1022 //
1023 // Revision 2.6 2009/11/23 23:44:32 ogrebeny
1024 // At Pibero's request, for the embedding infrastructure.
1025 //
1026 // Revision 2.5 2009/02/05 20:06:53 ogrebeny
1027 // Changed StEEmcDbMaker -> StEEmcDb
1028 //
1029 // Revision 2.4 2008/04/15 00:15:52 jwebb
1030 // Fixed bug in sector 12V.
1031 //
1032 // Revision 2.3 2008/04/11 14:37:17 jwebb
1033 // Added options to disable operation of individual slow simulaor subsystems.
1034 //
1035 // Revision 2.2 2007/11/28 16:17:30 jwebb
1036 // Added the following features:
1037 //
1038 // 1. User may specify the sampling fraction.
1039 //
1040 // 2. Tower and SMD gain spreads.
1041 //
1042 // Revision 2.1 2007/01/24 21:07:03 balewski
1043 // 1) no cout or printf, only new Logger
1044 // 2) EndcapMixer:
1045 // - no assert()
1046 // - locks out on first fatal error til the end of the job
1047 //
1048 // Revision 2.0 2007/01/13 00:03:04 jwebb
1049 // Upgrade of the slow simulator. The following changes have been made:
1050 //
1051 // 1. Towers will always be masked out when a "fail" bit is set in the database.
1052 // Previously this only happened if pedestals were being added, smeared.
1053 //
1054 // 2. Tower, preshower and postshower ADC values will be simulated using the
1055 // GEANT energy loss stored in StEmcHit and StMuEmcHit. Previously, ADC
1056 // values from the fast simulator were used and energy loss recovered
1057 // using gains, resulting in roundoff errors. Note that towers still use
1058 // the old path for MuDst-based analysis.
1059 //
1060 // 3. Tower simulation now accounts for the different light yields provided
1061 // by the brighter scintillator and two-fiber readout in the preshower
1062 // and postshower layers. Previously, only the difference in thickness
1063 // was accounted for by the GEANT simulation.
1064 //
1065 // Revision 1.5 2006/12/12 20:29:14 balewski
1066 // added hooks for Endcap embedding
1067 //
1068 // Revision 1.4 2006/08/07 18:50:11 balewski
1069 // added capabilty to run on StEvent, use se-method, see macros/ for example
1070 //
1071 // Revision 1.3 2005/09/23 17:23:01 balewski
1072 // now peds are added also to ADC of zero
1073 // fixed bug in ETOW ped smearing
1074 //
1075 // Revision 1.2 2005/09/23 01:30:10 jwebb
1076 // Tower peds now added if option is set.
1077 //
1078 // Revision 1.1 2004/12/15 17:02:56 balewski
1079 // try 2
1080 //
StMuEmcHit::getEnergy
float getEnergy() const
Return Hit energy.
Definition: StMuEmcHit.h:21
EEmcDbItem::name
char name[StEEmcNameLen]
ASCII name of the channel, see Readme.
Definition: EEmcDbItem.h:20
StEEmcSlowMaker::setAddPed
void setAddPed(Bool_t a=true)
Add pedestal offsets from DB.
Definition: StEEmcSlowMaker.h:255
StEEmcFastMaker::getPreshowerGain
static Float_t getPreshowerGain()
(adc=g*de ) fixed gain for pre/post shower
Definition: StEEmcFastMaker.cxx:333
StMaker::Make
virtual Int_t Make()
Definition: StMaker.cxx:898
StMuEmcHit::getAdc
int getAdc() const
Return ADC value.
Definition: StMuEmcHit.h:19
StEEmcFastMaker::getSmdGain
static Float_t getSmdGain()
(adc=g*de ) fixed gain for SMD
Definition: StEEmcFastMaker.cxx:326
StEEmcSlowMaker::setOverwrite
void setOverwrite(Bool_t o=true)
Overwrite the muDst values.
Definition: StEEmcSlowMaker.h:262
StEEmcSlowMaker::getMipdEdx
static Float_t getMipdEdx()
Return MIP dE/dx used for SMD, Pre, Post.
Definition: StEEmcSlowMaker.cxx:1001
StEEmcSlowMaker::setDropBad
void setDropBad(Bool_t d=true)
Drop bad channels marked as &quot;fail&quot; in DB.
Definition: StEEmcSlowMaker.h:259
StMuDst::muEmcCollection
static StMuEmcCollection * muEmcCollection()
returns pointer to current StMuEmcCollection
Definition: StMuDst.h:389
StEEmcSlowMaker::~StEEmcSlowMaker
virtual ~StEEmcSlowMaker()
Class destructor.
Definition: StEEmcSlowMaker.cxx:130
StMaker::GetName
virtual const char * GetName() const
special overload
Definition: StMaker.cxx:237
kStWarn
Definition: Stypes.h:42
StEEmcSlowMaker::Init
virtual Int_t Init()
Initialization.
Definition: StEEmcSlowMaker.cxx:157
StEEmcSlowMaker::setTowerGainSpread
void setTowerGainSpread(Float_t s, Float_t mean=1.0)
Defines a spread in the tower gains, generated gains will be between zero and mean + 1...
Definition: StEEmcSlowMaker.cxx:905
StEEmcSlowMaker::setSmearPed
void setSmearPed(Bool_t s=true)
Smear the pedestal with sigma from DB.
Definition: StEEmcSlowMaker.h:257
StEEmcSlowMaker
Slow simulator for EEMC.
Definition: StEEmcSlowMaker.h:137
StEEmcSlowMaker::setSource
void setSource(const Char_t *name)
Set the source of ADC. Can be &quot;MuDst&quot; (default) or &quot;StEvent&quot;.
Definition: StEEmcSlowMaker.cxx:603
kStErr
Definition: Stypes.h:44
StEEmcSlowMaker::setSmdGainSpread
void setSmdGainSpread(Float_t s, Int_t sector, Int_t uv, Int_t strip_index)
Defines a spread in the SMD gains, generated gains will be between zero and mean + 1...
Definition: StEEmcSlowMaker.cxx:936
kStOk
Definition: Stypes.h:41
StEEmcSlowMaker::Make
virtual Int_t Make()
Processes a single event.
Definition: StEEmcSlowMaker.cxx:238
Generated by   doxygen 1.8.5