StRoot: StEEmcSimulatorMaker/StEEmcSlowMaker.cxx Source File

00001 // *-- Author : Hal Spinka
00002 // 
00003 // $Id: StEEmcSlowMaker.cxx,v 2.11 2010/09/07 22:24:52 stevens4 Exp $
00004 
00005 #include <TFile.h>
00006 #include <TH2.h>
00007 #include <TRandom.h>
00008 #include <StMessMgr.h>
00009 
00010 #include "StEventTypes.h"
00011 #include "StMuDSTMaker/COMMON/StMuTypes.hh"
00012 
00013 #include "StEEmcUtil/database/EEmcDbItem.h"
00014 #include "StEEmcUtil/database/StEEmcDb.h"
00015 
00016 #include "StEEmcFastMaker.h"
00017 #include "StEEmcSlowMaker.h"
00018 
00019 ClassImp(StEEmcSlowMaker)
00020 
00021 //________________________________________________
00022 StEEmcSlowMaker::StEEmcSlowMaker(const Char_t *name, const Char_t*)
00023 : StMaker(name) {
00024    mMip2ene = getMipdEdx()*0.7; // This is the SMD thickness of 7 mm
00025                             // times the minimum ionizing energy loss of
00026                             // 1.998 MeV/cm from the PDG book
00027    mSig1pe = 0.85;          // from info from S. Vigdor on MAPMT test results
00028    //set different thicknesses for pre and post layers (found in geometry debug by Jason et. al.)
00029    mPmip2ene[0] = getMipdEdx()*0.475; // The pre-shower tiles are only 4.75 mm thick.
00030    mPmip2ene[1] = getMipdEdx()*0.475; // The pre-shower tiles are only 4.75 mm thick.
00031    mPmip2ene[2] = getMipdEdx()*0.5;   // The post-shower tiles are only 5 mm thick.
00032    mPmip2pe = 2.6*1.5;      // 2.6 mip/tower scint * 1.5 light yield
00033                             // in pre- and post-shower elements
00034    // loop to init  mMip2pe[] - this will eventually need to be
00035    // replaced by a table of measured values
00036    for (Int_t i = 0;i < MaxSmdStrips;i++) {
00037      // This is currently a trapezoidal-parameterization of
00038      // the light curve measurements
00039      mMip2pe[i] = avgNumPePerMip(i);
00040    }
00041 
00042   mEeDb=0;
00043   mNInpEve=0; 
00044   memset(mHist,0,sizeof(mHist));
00045 
00047   mEnableTower=true;
00048   mEnableSMD=true;
00049   mEnablePrePost=true;
00050 
00052   mAddPed   = true;                                  
00054   mSmearPed = true;                                  
00056   mDropBad  = false;
00058   mOverwrite = true;
00059 
00060   mIsEmbeddingMode = false;
00062   mSource = kMuDst;
00063 
00065   mTruncatePedSmear = 3;
00066   
00069   setDoLightYield(true);
00070 
00077   setRelativeLightYield( (4.0/4.75) * 1.68, (4.0/4.75) * 1.68, (4.0/5.0)* 0.94);
00078 
00080   mSamplingFraction = StEEmcFastMaker::getSamplingFraction();
00081   mSamplingFractionUser = mSamplingFraction;
00082 
00084   for (Int_t i = 0;i < kEEmcNumEtas;i++) {
00085       mTowerGains[i] = StEEmcFastMaker::getTowerGains()[i] * mSamplingFraction;
00086   }
00087   mPrepostGains = StEEmcFastMaker::getPreshowerGain();
00088   mSmdGains     = StEEmcFastMaker::getSmdGain();
00089   mMaxAdc       = StEEmcFastMaker::getMaxAdc();
00090   
00091   if ( (mPrepostGains-23000.0)>1. ||
00092        (mSmdGains-23000.0)>1.     ||
00093        (mTowerGains[0]-18.9654)>0.01 )
00094     {
00095       LOG_WARN << "--------------------------------------------" << endm;
00096       LOG_WARN << "You changed the gains in the fast simulator." << endm;
00097       LOG_WARN << "You are on your own." << endm;
00098       LOG_WARN << "                   -- the eemc software team" << endm;
00099       LOG_WARN << endm;
00100       LOG_WARN << "mSmdGains=" << mSmdGains << endm;
00101       LOG_WARN << "mPrepostGains="<<mPrepostGains<<endm;
00102       for ( Int_t ii=0;ii<12;ii++ )
00103         LOG_WARN << "mTowerGains[etabin="<<ii<<"]="<<mTowerGains[ii]<<endm;
00104       //assert(2+2==5);
00105     }
00106 
00107   // initialize tower gain factors to 1
00108   for (Int_t sec = 0;sec < kEEmcNumSectors;sec++) {
00109     for (Int_t sub = 0;sub < kEEmcNumSubSectors;sub++) {
00110       for (Int_t eta = 0;eta < kEEmcNumEtas;eta++) {
00111           mTowerGainFact[sec][sub][eta] = 1.0;
00112       }
00113     }
00114   }
00115 
00116   // initialize SMD gain factors to 1
00117   for (Int_t sec = 0;sec < kEEmcNumSectors;sec++) {
00118     for (Int_t uv = 0;uv < kEEmcNumSmdUVs;uv++) {
00119       for (Int_t strip = 0;strip < kEEmcNumStrips;strip++) {
00120           mSmdGainFact[sec][uv][strip] = 1.0;
00121       }
00122     }
00123   }
00124 
00125   // flag set if used in BFC to change defaults in Init()
00126   mIsBFC = GetParentChain()->InheritsFrom("StBFChain");         
00127 }
00128 
00129 //________________________________________________
00130 StEEmcSlowMaker::~StEEmcSlowMaker() {
00131 }
00132 
00133 //________________________________________________
00134 Float_t StEEmcSlowMaker::avgNumPePerMip(Int_t stripID) {                                                                                                           
00135 //                                                                                                                                                          
00136 // A parameterization of the average number of photoelectrons                                                                                               
00137 // per mip for a given SMD strip.  See elog 457.                                                                                                            
00138 //                                                                                                                                                          
00139   Float_t ya=0,yb=0,xa=1,xb=2;                                                                                                                                 
00140   if ( stripID<1) {                                                                                                                                          
00141     ;                                                                                                                                                        
00142   } else if (stripID<20) {                                                                                                                                   
00143     xa=1; ya=2.;                                                                                                                                             
00144     xb=20; yb=4.;                                                                                                                                            
00145   } else if (stripID<250) {                                                                                                                                  
00146     xa=20; ya=4.;                                                                                                                                            
00147     xb=250; yb=6.;                                                                                                                                           
00148   } else {                                                                                                                                                   
00149     xa=250; ya=6.;                                                                                                                                           
00150     xb=290; yb=9.;                                                                                                                                           
00151   }                                                                                                                                                          
00152   const Float_t y = ya + (yb - ya) / (xb - xa) * (stripID - xa);
00153   return y;
00154 }                                                                                                                                                            
00155  
00156 //________________________________________________
00157 Int_t StEEmcSlowMaker::Init() {
00158   LOG_INFO << "mIsEmbeddingMode=" << mIsEmbeddingMode << ", mIsBFC=" << mIsBFC << endm;          
00159 
00160   if (mIsEmbeddingMode) {
00161     setDropBad(1);   // force bad channels to be dropped in db
00162     setAddPed(0);    // disable pedestal addition
00163     setSmearPed(0);  // disable pedestal smearing
00164     setOverwrite(1); // overwrites StEvent values
00165     setSource("StEvent");    
00166   } else if (mIsBFC) { //set defaults for BFC with no embedding 
00167     setAddPed(1);    // add pedestals                 
00168     setSmearPed(1);  // smear pedestals                 
00169     setDropBad(1);   // force bad channels to be dropped in db
00170     setOverwrite(1); // overwrite ADC values
00171     setSource("StEvent");             
00172   } else { //set defaults for running in analysis chain 
00173     if (mSource == kMuDst) disableTower(); //don't change tower in analysis chain since only ADC is stored in MuDst  
00174   }                                   
00175     
00176   //print out full configuration in log file  
00177   LOG_INFO << "mAddPed=" << mAddPed << ", mSmearPed=" << mSmearPed << ", mDropBad=" << mDropBad << ", mOverwrite=" << mOverwrite 
00178     << ", mEnableTower=" << mEnableTower << ", mEnableSMD=" << mEnableSMD << ", mEnablePrePost=" << mEnablePrePost << endm;
00179 
00180   mEeDb = (StEEmcDb*)GetDataSet("StEEmcDb");
00181   if (!mEeDb) {
00182     LOG_ERROR << "Cannot find EEMC database StEEmcDb" << endm;
00183   }
00184   InitHisto();
00185   if ( mSmearPed && !mAddPed) {
00186     LOG_WARN << "detected mSmearPed=true && mAddPed=false\nwill not work due to ETOW hits storage container in muDst not accepting negative ADC values." << endm; 
00187     //assert(0);
00188   }
00189   if ( mSmearPed ) {
00190     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; 
00191   }
00192   return StMaker::Init();
00193 }
00194 
00195 //________________________________________________
00196 void StEEmcSlowMaker::InitHisto() {
00197   Char_t tt1[100], tt2[500];
00198 
00199   sprintf(tt1,"mm");
00200   sprintf(tt2,"freq vs. sector ID; sector ID");
00201   mHist[12]=new TH1F(tt1,tt2,20,-0.5,19.5);
00202   
00203   sprintf(tt1,"PreADC");
00204   sprintf(tt2,"Revised ADC for Pre/Post Scint.");
00205   mHist[0]=new TH1F(tt1,tt2,100,0,200.);
00206   
00207   sprintf(tt1,"Pre2ADC");
00208   sprintf(tt2,"Initial ADC for Pre/Post Scint.");
00209   mHist[1]=new TH1F(tt1,tt2,100,0,200.);
00210   
00211   sprintf(tt1,"SMDADC");
00212   sprintf(tt2,"Revised ADC for SMD Strips");
00213   mHist[2]=new TH1F(tt1,tt2,100,0,200.);
00214 
00215   sprintf(tt1,"SMD2ADC");
00216   sprintf(tt2,"Initial ADC for SMD Strips");
00217   mHist[3]=new TH1F(tt1,tt2,100,0,200.);
00218 
00219   mHist[4]=new TH1F("hADCtow","Initial ADC for towers",512+32,-32.,512.);
00220   mHist[5]=new TH1F("hADCtow2","Finial ADC fot towers",512+32,-32.,512.);
00221 
00222   //..................
00223   sprintf(tt1,"adc");
00224   sprintf(tt2," Adc vs strip ID; strip ID; ADC ");
00225   mHist[20]=new TH2F(tt1,tt2,30,0,300,100,0,200.);
00226 
00227   sprintf(tt1,"ADCvsstr");
00228   sprintf(tt2," Adc vs, strip ID; strip ID; ADC ");
00229   mHist[19]=new TH2F(tt1,tt2,60,0,300,100,0,200.);
00230 
00231   // add histos to the list (if provided)
00232   for (Int_t i = 0; i < maxHist; ++i) {
00233       if (mHist[i]) this->AddHist(mHist[i]);
00234   }
00235 }
00236 
00237 //________________________________________________
00238 Int_t StEEmcSlowMaker::Make() {
00239   Int_t result = StMaker::Make();
00240   mNInpEve++;
00241   LOG_DEBUG << "iEve " << mNInpEve << ", mSource = " << mSource << endm;
00242   
00243   switch (mSource) {
00244   case kMuDst:
00246     {
00247       if (!GetInputDS("MuDst")) {
00248         LOG_DEBUG<<"::Make() MuDst missing"<<endm;
00249         return kStWarn;
00250       }
00251 
00252       StMuEmcCollection* emc = StMuDst::muEmcCollection();
00253       if (!emc) {
00254         LOG_DEBUG<<"::Make() StMuEmcCollection missing"<<endm;
00255         return kStWarn;
00256       }
00257 
00259       if (mEnableTower && (result == kStOk)) result = MakeTower(emc);
00260 
00262       if (mEnablePrePost && (result == kStOk)) result = MakePrePost(emc);
00263 
00265       if (mEnableSMD && (result == kStOk)) result = MakeSMD(emc);
00266     }
00267     break;
00268 
00269   case kStEvent:
00271     {
00272 
00273       StEmcCollection *emc =0; 
00274       if(mIsEmbeddingMode) {
00275         StEEmcFastMaker *fast = (StEEmcFastMaker*)GetMakerInheritsFrom("StEEmcFastMaker");
00276         if(fast==0) {
00277           LOG_WARN << GetName() << "::Make()  no EEmcFastSim in the chain, ignore Endcap"<< endm;
00278           return kStWarn;
00279         }
00280         emc = fast->GetLocalEmcCollection();
00281       } else { // it is not Embedding mode
00282         
00283         StEvent* event = (StEvent*)GetInputDS("StEvent");
00284         if (!event) {
00285           LOG_WARN << GetName() << "::Make()  no StEvent"<< endm;
00286           return kStWarn;
00287         }
00288         emc = event->emcCollection();
00289       }
00290       // now emc-collection should be accessible one way or another
00291       
00292       if (!emc) {
00293         LOG_WARN << GetName() << "::Make()  no emcCollection()" << endm;
00294         return kStWarn;
00295       }
00296 
00298       if (mEnableTower && (result == kStOk)) result = MakeTower(emc);
00299 
00301       if (mEnablePrePost && (result == kStOk)) result = MakePrePost(emc);
00302 
00304       if (mEnableSMD && (result == kStOk)) result = MakeSMD(emc);
00305     }
00306     break;
00307 
00308   default:
00309     LOG_ERROR<< "Unknown source type " << mSource << " for this event" << endm;
00310     break;
00311   }
00312   return result;
00313 }
00314 
00315 //________________________________________________
00316 //
00317 // Some notes on differences between using StEvent and StMuDst
00318 //
00319 //    Roundoff errors
00320 // 
00321 //    Since the MuDst only stores ADC values for the towers, it will be 
00322 //    subject to roundoff errors.  This is most noticable for towers with
00323 //    small ADC response.  Keep in mind that (int)2.7 rounds down to 2.
00324 //
00325 //    To partially offset this, we will always assume a starting ADC value
00326 //    of ADC+0.5 for the MuDst.
00327 //
00328 
00329 Int_t StEEmcSlowMaker::MakeTower(StMuEmcCollection *emc) {
00330 
00331   /**************************************************************
00332    *
00333    * Performs pedestal addition and smearing, computes ADC based
00334    * on GEANT energy deposit and database gains, and corrects
00335    * tower ADC for the relative brightness of the pre- and post-
00336    * shower layers.
00337    *
00338    **************************************************************
00339    */
00340 
00342   Float_t prepost_adc[kEEmcNumSectors][kEEmcNumSubSectors][kEEmcNumEtas];
00343   memset( prepost_adc, 0, sizeof(prepost_adc) );
00344  
00348   if (mDoLightYield) for (Int_t i = 0;i < emc->getNEndcapPrsHits();i++) {
00350       Int_t sec,sub,eta,pre;
00351       StMuEmcHit *hit=emc->getEndcapPrsHit(i,sec,sub,eta,pre);
00352       // range check on returned values
00353       if (!( sec >= 1 && sec <= 12) || !(sub >= 1 && sub <= 5) || !(eta >= 1 && eta <= 12) || !(pre >= 1 && pre <= 3)) {
00354         LOG_ERROR << "Indexing errors detected: EPRS hit " << i << ", sec = " << sec << ", sub = " << sub << ", eta = " << eta << ", pre = " << pre << endm;
00355         setZeroAdc(emc);
00356         return kStErr;
00357       }
00360       const Float_t edeposit = hit ? hit->getEnergy() : 0;
00361 
00362       // divide by sampling fraction, multiply by _tower_ gain
00363       // and brightness factor, and add to the adc sum for 
00364       // this tower.  Since GEANT already accounts for the 
00365       // difference in thickness between pre/post and normal
00366       // layers, the factor of 0.8 is introduced to prevent us
00367       // from double correcting.
00368       const EEmcDbItem *tower = mEeDb ? mEeDb->getTile(sec,sub-1+'A', eta, 'T') : 0;
00369       if (!tower) {
00370         LOG_ERROR << "Cannot find DB entry for ETOW: sec = " << sec << ", sub = " << sub << ", eta = " << eta << endm;
00371         continue;
00372       }
00373       Float_t myadc = edeposit / mSamplingFractionUser * tower->gain;
00374       myadc *= ( mRelativeLightYield[pre-1] - 1.0 );
00375 
00376       prepost_adc[sec-1][sub-1][eta-1] += myadc;
00377   }
00378 
00382   for (Int_t i = 0;i < emc->getNEndcapTowerADC();i++) {
00383       Int_t sec,sub,eta,adc;
00385       emc->getEndcapTowerADC(i,adc,sec,sub,eta);
00386       // range check on returned values
00387       if (!(sec >= 1 && sec <= 12) || !(sub >= 1 && sub <= 5) || !(eta >= 1 && eta <= 12)) {
00388         LOG_ERROR << "Indexing errors detected: ETOW hit " << i << ", sec = " << sec << ", sub = " << sub << ", eta = " << eta << endm;
00389         setZeroAdc(emc);
00390         return kStErr;
00391       }
00392       
00393       const Int_t old = adc;
00394       if (mHist[4]) mHist[4]->Fill( old );
00395 
00397       const EEmcDbItem *tower = mEeDb ? mEeDb->getTile(sec,sub-1+'A', eta, 'T') : 0;
00398       if (!tower) {
00399         LOG_ERROR << "Cannot find DB entry for ETOW: sec = " << sec << ", sub = " << sub << ", eta = " << eta << endm;
00400         continue;
00401       }
00402 
00403       if (mDropBad && (tower->fail || !checkDBped(tower))) {
00404         //zero out values if tower has fail bit set or ped is bad        
00405         adc = 0;
00406       } else {
00408         Float_t myadc=(Float_t)adc + 0.5;
00411         myadc *= tower->gain / mTowerGains[eta-1] * mSamplingFraction / mSamplingFractionUser;
00413         myadc += prepost_adc[sec-1][sub-1][eta-1];
00415         myadc *= ( mTowerGainFact[sec-1][sub-1][eta-1] );
00416         Float_t ped = tower->ped;
00418         if (mAddPed) {
00419             if ( mSmearPed ) ped += getPedSmear( tower->sigPed );   
00420             myadc += ped;
00421         }
00423         adc = (Int_t)myadc;
00424         if (mHist[5]) mHist[5]->Fill(adc);
00425         if (mDropBad && (tower->gain <= 0) && mAddPed) {
00426             //still add ped to tower if gain is bad since trigger emulator expects peds
00427             adc = ped;
00428         }
00430         if (adc < 0) adc = 0;
00431         if (adc > mMaxAdc) adc = mMaxAdc;
00432       }
00433       if (mOverwrite) {
00434         if (old) {
00435             LOG_DEBUG << "overwriting tower=" << tower->name << " old adc=" << old << " new adc=" << adc << endm;
00436         }
00442         emc->setTowerADC( i+1, adc, eemc );
00443       }
00444   }
00445   return kStOk;
00446 }
00447 
00448 //________________________________________________
00449 Int_t StEEmcSlowMaker::MakePrePost(StMuEmcCollection *emc) {
00450   for (Int_t i = 0;i < emc->getNEndcapPrsHits();i++) {
00451     Int_t pre,sec,eta,sub;
00453     StMuEmcHit *hit = emc->getEndcapPrsHit(i,sec,sub,eta,pre);
00454     if (!hit) continue;
00455     
00456     // range check on returned values
00457     if (!( sec >= 1 && sec <= 12) || !(sub >= 1 && sub <= 5) || !(eta >= 1 && eta <= 12) || !(pre >= 1 && pre <= 3)) {
00458       LOG_ERROR << "Indexing errors detected: EPRS hit " << i << ", sec = " << sec << ", sub = " << sub << ", eta = " << eta << ", pre = " << pre << endm;
00459       setZeroAdc(emc);
00460       return kStErr;
00461     }
00462 
00464     if (mEeDb && (sec < mEeDb->getFirstSector() || sec > mEeDb->getLastSector())) continue;
00465      
00467     const EEmcDbItem *x = mEeDb ? mEeDb->getTile(sec,sub-1+'A', eta, pre-1+'P') : 0; 
00468     if (!x) {
00469         LOG_ERROR << "Cannot find DB entry for EPRS: sec = " << sec << ", sub = " << sub << ", eta = " << eta << ", pre = " << pre << endm;
00470         continue;
00471     }
00472 
00473     Int_t NUadc;
00474     // Zero out values if tower is marked as bad
00475     if (mDropBad && (x->fail || !checkDBped(x) || (x->gain <= 0))) {
00476         NUadc=0;
00477     } else {
00478         const Float_t adc = hit->getAdc();
00479         const Float_t energy = hit->getEnergy();
00480 
00481         Float_t newadc = energy * x->gain;
00482         if (adc > 0) {
00483             if (mHist[1]) mHist[1]->Fill(adc);
00484       
00487             const Float_t mipval = energy / mPmip2ene[pre-1];
00488             const Float_t avgpe  = mipval * mPmip2pe;
00489             const Float_t Npe    = gRandom->Poisson(avgpe);
00490       
00491             if (Npe > 0) {
00493                 const Float_t sigmape   = sqrt(Npe) * mSig1pe;
00494                 Float_t smearedpe = gRandom->Gaus(Npe,sigmape);
00495                 if (smearedpe < 0) smearedpe = 0;
00497                 newadc = smearedpe * mMip2ene*(x->gain) / mPmip2pe;
00498             }
00499         } // non-zero ADC on input
00500 
00502         Float_t ped = (mAddPed) ? x->ped : 0;
00504         if ( mSmearPed ) ped += getPedSmear( x->sigPed ); 
00505        
00507         NUadc = (Int_t)(newadc + 0.5 + ped );
00508     
00510         if (NUadc < 0) NUadc = 0;
00511         if (NUadc > mMaxAdc) NUadc = mMaxAdc;
00512         if (adc > 0 && mHist[0]) mHist[0]->Fill(NUadc); //??
00513     }
00514 
00518     if (mOverwrite) hit->setAdc(NUadc);
00519   }
00520   return kStOk;
00521 }
00522 
00523 //________________________________________________
00524 Int_t StEEmcSlowMaker::MakeSMD(StMuEmcCollection *emc) {
00525   Int_t iuv = 0;
00526   for (Char_t uv = 'U';uv <= 'V';uv++) {
00527     iuv++;
00528     Int_t sec,strip;
00529     for (Int_t i = 0;i < emc->getNEndcapSmdHits(uv);i++) {
00530       StMuEmcHit *hit=emc->getEndcapSmdHit(uv,i,sec,strip);
00531       // range check on returned values
00532       if (!(strip-1 >= 0) || !(strip <= 288) || !(iuv-1 >= 0 && iuv-1 < 2) || !(sec > 0 && sec <= MaxSectors)) {
00533         LOG_ERROR << "Bad index for ESMD: sec = " << sec << ", iuv = " << iuv << ", strip = " << strip << endm;
00534         setZeroAdc(emc); 
00535         return kStErr;
00536       }
00537       
00538       // tmp, for fasted analysis use only hits from sectors init in DB
00539       if (mEeDb && (sec < mEeDb->getFirstSector() || sec > mEeDb->getLastSector())) continue;
00540       const EEmcDbItem *x = mEeDb ? mEeDb->getByStrip(sec,uv,strip) : 0;
00541       if (!x) {
00542         LOG_ERROR << "Cannot find DB entry for ESMD: sec = " << sec << ", uv = " << uv << ", strip = " << strip << endm;
00543         continue;
00544       }
00545 
00546       Int_t NUadc;
00547       // Zero out values if tower is marked as bad
00548       if (mDropBad && (x->fail || !checkDBped(x) || (x->gain <= 0))) {
00549           NUadc=0;
00550       } else {
00551         const Float_t adc = hit->getAdc();
00552         const Float_t energy = hit->getEnergy();
00553         Float_t newadc = energy * x->gain;
00554       
00555         if (adc > 0) {
00556             if (mHist[12]) mHist[12]->Fill(x->sec);
00557             if (mHist[20]) mHist[20]->Fill(x->strip,adc); 
00558             if (mHist[3]) mHist[3]->Fill(adc);
00559 
00560             // Addition of Poisson fluctuations and 
00561             // Gaussian 1 p.e. resolution
00562             const Float_t mipval = energy / mMip2ene;
00563             const Float_t avgpe  = mipval * mMip2pe[strip];
00564             const Float_t Npe    = gRandom->Poisson(avgpe);
00565             if (Npe > 0) {
00567                 const Float_t sigmape   = sqrt(Npe) * mSig1pe;
00568                 Float_t smearedpe = gRandom->Gaus(Npe, sigmape);
00569                 if (smearedpe < 0) smearedpe = 0;
00571                 newadc = smearedpe * mMip2ene * (x->gain) / mMip2pe[strip];
00572             }
00573         }
00574       
00575         // Add in factor from gain smearing
00576         newadc *= mSmdGainFact[sec-1][iuv-1][strip-1];
00577 
00578         // Lookup pedestal in database (possibly zero)
00579         Float_t ped = (mAddPed) ? x->ped : 0;
00581         if ( mSmearPed ) ped += getPedSmear( x->sigPed );
00582             
00584         NUadc = (Int_t)(newadc + 0.5 + ped );
00585       
00587         if (NUadc < 0) NUadc = 0;
00588         if (NUadc > mMaxAdc) NUadc = mMaxAdc;
00589         if (adc > 0 && mHist[19]) mHist[19]->Fill(x->strip,NUadc);
00590         if (adc > 0 && mHist[2]) mHist[2]->Fill(NUadc);
00591       }
00592       
00596       if (mOverwrite) hit->setAdc(NUadc);
00597     }// loop over 1 plane
00598   } // loop over U,V
00599   return kStOk;
00600 }
00601 
00602 //________________________________________________
00603 void StEEmcSlowMaker::setSource(const Char_t* name) {
00604   if (strcmp(name, "MuDst") == 0) {
00605     mSource = kMuDst;
00606   } else if (strcmp(name, "StEvent") == 0) {
00607     mSource = kStEvent;
00608   } else {
00609     LOG_WARN<<"::setSource()"<<"Source must be \"MuDst\" or \"StEvent\""<<endm;
00610   }
00611 }
00612 
00613 //________________________________________________
00614 Int_t StEEmcSlowMaker::MakeTower(StEmcCollection* emc) {
00615   StEmcDetector *det = emc->detector(kEndcapEmcTowerId);
00616   if (!det) {
00617     LOG_DEBUG<<"No kEndcapEmcTowerId"<<endm;
00618     return kStOk;
00619   } 
00620   StEmcDetector *prepost = emc->detector(kEndcapEmcPreShowerId);
00621   if (!prepost) {
00622     LOG_WARN<<"No kEndcapEmcPreShowerId, towers are fast simu only."<<endm;
00623     return kStOk; 
00624   } 
00625 
00627   Float_t prepost_adc[kEEmcNumSectors][kEEmcNumSubSectors][kEEmcNumEtas];
00628   memset( prepost_adc, 0, sizeof(prepost_adc) );
00629 
00630   for (UInt_t sec = 1;sec <= det->numberOfModules();sec++) {
00631       StSPtrVecEmcRawHit &prepost_hits = prepost->module(sec)->hits();
00632       if (mDoLightYield) for (UInt_t ihit = 0;ihit < prepost_hits.size();ihit++) {
00633           StEmcRawHit *hit = prepost_hits[ihit];
00634           const UInt_t sub = (hit->sub()-1)%5+1;
00635           const UInt_t eta = hit->eta();
00636           const UInt_t pre = (hit->sub()-1)/5+1;
00637           // range check on returned values
00638           if (!(sec >= 1 && sec <= 12) || !(sub >= 1 && sub <= 5) || !(eta >= 1 && eta <= 12) || !(pre >= 1 && pre <= 3)) {
00639             LOG_ERROR << "Indexing errors detected for EPRS: sec = " << sec << ", sub = " << sub << ", eta = " << eta << ", pre = " << pre << endm;
00640             setZeroAdc(emc); 
00641             return kStErr;
00642           }
00643                   
00644           // energy deposit in GeV
00645           const Float_t edeposit = hit->energy();
00646           // divide by sampling fraction, multiply by _tower_ gain
00647           // and brightness factor, and add to the adc sum for 
00648           // this tower.  Since GEANT already accounts for the 
00649           // difference in thickness between pre/post and normal
00650           // layers, the factor of 0.8 is introduced to prevent us
00651           // from double correcting.
00652           const EEmcDbItem *tower = mEeDb ? mEeDb->getTile(sec,sub-1+'A', eta, 'T') : 0;
00653           if (!tower) {
00654             LOG_ERROR << "Cannot find DB entry for ETOW: sec = " << sec << ", sub = " << sub << ", eta = " << eta << endm;
00655             continue;
00656           }
00657           Float_t myadc = edeposit / mSamplingFraction * tower->gain;
00658           myadc *= ( mRelativeLightYield[pre-1] - 1.0 );
00659           prepost_adc[sec-1][sub-1][eta-1] += myadc;
00660       }
00661       Int_t nTchange = 0; // counts # of towers with ADC (int) values changed
00662       StSPtrVecEmcRawHit &tower_hits = det->module(sec)->hits();
00663       for (UInt_t ihit = 0;ihit < tower_hits.size();ihit++) {
00664           StEmcRawHit *hit = tower_hits[ihit];
00665           const UInt_t sub = hit->sub();
00666           const UInt_t eta = hit->eta();
00667           // range check on returned values
00668           if (!(sec >= 1 && sec <= 12) || !(sub >= 1 && sub <= 5) || !(eta >= 1 && eta <= 12)) {
00669             LOG_ERROR << "Indexing errors detected for ETOW: sec = " << sec << ", sub = " << sub << ", eta = " << eta << endm;
00670             setZeroAdc(emc); 
00671             return kStErr;
00672           }
00673 
00674           // get DB entry for this tower
00675           const EEmcDbItem *tower = mEeDb ? mEeDb->getTile(sec,sub-1+'A', eta, 'T') : 0;
00676           if (!tower) {
00677             LOG_ERROR << "Cannot find DB entry for ETOW: sec = " << sec << ", sub = " << sub << ", eta = " << eta << endm;
00678             continue;
00679           }
00680           Int_t adc;
00681           if (mDropBad && (tower->fail || !checkDBped(tower))) {
00682             //zero out values if tower has fail bit set or ped is bad
00683             adc = 0;
00684           } else {
00685             if (mHist[4]) mHist[4]->Fill( hit->adc() );
00686             // energy deposit in GeV
00687             const Float_t edeposit = hit->energy();
00688             // compute adc using gains stored in the databse
00689             Float_t myadc = edeposit / mSamplingFraction * tower->gain;
00690             // add in brightness correction for the pre and postshower layers
00691             myadc += prepost_adc[sec-1][sub-1][eta-1];
00693             myadc *= ( mTowerGainFact[sec-1][sub-1][eta-1] );
00695             Float_t ped = tower->ped;
00696             if (mAddPed) {
00697                 if ( mSmearPed ) ped += getPedSmear( tower->sigPed );
00698                 myadc += ped;
00699             }
00700             adc = (Int_t)myadc;
00701             // Zero out values if tower is marked as bad
00702             if (mDropBad && (tower->gain <= 0) && mAddPed) {
00703                 //still add ped to tower if gain is bad since trigger emulator expects peds
00704                 adc = ped;
00705             }
00707             if (adc < 0) adc = 0;
00708             if (adc > mMaxAdc) adc = mMaxAdc;
00709           }
00710           if (mHist[5]) mHist[5]->Fill(adc);
00711           if (mOverwrite) {
00712             if ((Int_t)hit->adc() != adc) nTchange++;
00713             LOG_DEBUG <<"overwriting tower=" << tower->name << " old adc=" << hit->adc() << " new adc=" << adc << endm;
00714             hit->setAdc( adc );
00715           }
00716       }
00717       LOG_DEBUG << " Etow changed " << nTchange << " ADC values for sector=" << sec << endm;
00718   }// loop over sectors
00719   return kStOk;
00720 }
00721 
00722 //________________________________________________
00723 Int_t StEEmcSlowMaker::MakePrePost(StEmcCollection* emc) {
00724   StEmcDetector* det = emc->detector(kEndcapEmcPreShowerId);
00725   if (!det) {
00726     LOG_WARN << "No kEndcapEmcPreShowerId" << endm;
00727     return kStOk;
00728   }
00729 
00730   for (UInt_t sector =1; sector <= det->numberOfModules(); ++sector) {
00731     StSPtrVecEmcRawHit& hits = det->module(sector)->hits();
00732     for(UInt_t i = 0; i < hits.size(); ++i) {
00733       StEmcRawHit* hit = hits[i];
00734       Char_t sub = 'A'+(hit->sub()-1)%5;
00735       
00736       // Layer: 'P'=preshower1, 'Q'=preshower2, 'R'=postshower
00737       Char_t layer = 'P'+(hit->sub()-1)/5;
00738       Int_t layerP = (hit->sub()-1)/5;
00739 
00740       const UInt_t isub = (hit->sub()-1)%5+1;
00741       const UInt_t ieta = hit->eta();
00742       const UInt_t ipre = (hit->sub()-1)/5+1;
00743       // range check on returned values
00744       if (!(sector >= 1 && sector <= 12) || !(isub >= 1 && isub <= 5) || !(ieta >= 1 && ieta <= 12) || !(ipre >= 1 && ipre <= 3)) {
00745         LOG_ERROR << "Indexing errors detected for EPRS: sec = " << sector << ", sub = " << isub << ", eta = " << ieta << ", pre = " << ipre << endm;
00746         setZeroAdc(emc); 
00747         return kStErr;
00748       }
00749       
00750       // Database ranges: sector=1-12, sub=A-E, eta=1-12, type=T,P-R; Slow method
00751       const EEmcDbItem* x = mEeDb ? mEeDb->getTile(sector, sub, hit->eta(), layer) : 0;
00752       if (!x) {
00753         LOG_ERROR << "Cannot find DB entry for EPRS: sec = " << sector << ", sub = " << sub << ", eta = " << hit->eta() << ", pre = " << layer << endm;
00754         continue;
00755       }
00756 
00757       Int_t NUadc;
00758       // Zero out values if tower is marked as bad
00759       if (mDropBad && (x->fail || !checkDBped(x) || (x->gain <= 0))) {
00760           NUadc=0;
00761       } else {
00762         Float_t adc = hit->adc();
00763         Float_t energy = hit->energy();
00764         Float_t newadc = energy * x->gain;
00765         if(adc > 0) {
00766             if (mHist[1]) mHist[1]->Fill(adc);
00767 
00770             const Float_t mipval = energy / mPmip2ene[layerP];
00771             const Float_t avgpe  = mipval * mPmip2pe;
00772             const Float_t Npe    = gRandom->Poisson(avgpe);
00773 
00774             if (Npe > 0) {
00776                 const Float_t sigmape   = sqrt(Npe) * mSig1pe;
00777                 Float_t smearedpe = gRandom->Gaus(Npe, sigmape);
00778                 if (smearedpe < 0) smearedpe = 0;
00780                 newadc = smearedpe * mMip2ene * x->gain / mPmip2pe;
00781             }
00782         } // non-zero ADC on input
00783 
00785         Float_t ped = mAddPed ? x->ped : 0;
00786 
00788         if ( mSmearPed ) ped += getPedSmear( x->sigPed );
00789       
00791         NUadc = Int_t(newadc + 0.5 + ped);
00792 
00794         if ( NUadc < 0 ) NUadc = 0;
00795         if ( NUadc > mMaxAdc ) NUadc = mMaxAdc;
00796         if (adc > 0 && mHist[0]) mHist[0]->Fill(NUadc); //??
00797       }
00798 
00802       if (mOverwrite) hit->setAdc(NUadc);
00803     }
00804   }
00805   return kStOk;
00806 }
00807 
00808 //________________________________________________
00809 Int_t StEEmcSlowMaker::MakeSMD(StEmcCollection* emc) {
00810   Int_t iuv = 0;
00811   for (Char_t plane = 'U'; plane <= 'V'; ++plane) {
00812     iuv++;
00813     StEmcDetector* det = 0;
00814     switch (plane) {
00815     case 'U':
00816       det = emc->detector(kEndcapSmdUStripId);
00817       if (!det) {
00818         LOG_DEBUG << "No kEndcapSmdUStripId" << endm;
00819         continue;
00820       }
00821       break;
00822     case 'V':
00823       det = emc->detector(kEndcapSmdVStripId);
00824       if (!det) {
00825         LOG_DEBUG << "No kEndcapSmdVStripId" << endm;
00826         continue;
00827       }
00828       break;
00829     }
00830 
00831     if (!det) continue;
00832 
00833     for (UInt_t sector = 1; sector <= det->numberOfModules();++sector) {
00834       StSPtrVecEmcRawHit& hits = det->module(sector)->hits();
00835       for (UInt_t i = 0; i < hits.size(); ++i) {
00836         StEmcRawHit* hit = hits[i];
00837         Int_t strip = hit->eta();
00838         
00839         // range check on returned values
00840         if (!(strip-1 >= 0) || !(strip <= 288) || !(iuv-1 >= 0 && iuv-1 < 2) || !(sector > 0 && sector <= MaxSectors)) {
00841           LOG_ERROR << "Bad index for ESMD: sec = " << sector << ", iuv = " << iuv << ", strip = " << strip << endm;
00842           setZeroAdc(emc); 
00843           return kStErr;
00844         }
00845 
00846         // Database ranges: sector=1-12, plane=U-V, strip=1-288
00847         const EEmcDbItem* x = mEeDb ? mEeDb->getByStrip(sector, plane, strip) : 0;
00848         if (!x) {
00849           LOG_ERROR << "Cannot find DB entry for ESMD: sec = " << sector << ", uv = " << iuv-1+'U' << ", strip = " << strip << endm;
00850           continue;
00851         }
00852 
00853         Int_t NUadc;
00854         // Zero out values if tower is marked as bad
00855         if (mDropBad && (x->fail || !checkDBped(x) || (x->gain <= 0))) {
00856             NUadc = 0;
00857         } else {
00858             const Float_t adc    = hit->adc();
00859             const Float_t energy = hit->energy();
00860             Float_t newadc = energy * x->gain;
00861             if (adc > 0) {
00862                 if (mHist[12]) mHist[12]->Fill(x->sec);
00863                 if (mHist[2]) mHist[20]->Fill(x->strip, adc); 
00864                 if (mHist[3]) mHist[3]->Fill(adc);
00865                 // Addition of Poisson fluctuations and 
00866                 // Gaussian 1 p.e. resolution
00867                 const Float_t mipval = energy / mMip2ene;
00868                 const Float_t avgpe  = mipval * mMip2pe[strip];
00869                 const Float_t Npe    = gRandom->Poisson(avgpe);
00870                 if (Npe > 0) {
00872                     const Float_t sigmape   = sqrt(Npe) * mSig1pe;
00873                     Float_t smearedpe = gRandom->Gaus(Npe, sigmape);
00874                     if (smearedpe < 0) smearedpe = 0;
00876                     newadc = smearedpe * mMip2ene * x->gain / mMip2pe[strip];
00877                 }
00878             }
00879             // Add in factor from gain smearing
00880             newadc *= mSmdGainFact[sector-1][iuv-1][strip-1];
00882             Float_t ped = mAddPed ? x->ped : 0;
00884             if ( mSmearPed ) ped += getPedSmear( x->sigPed );
00886             NUadc = Int_t(newadc + 0.5 + ped);
00888             if ( NUadc < 0 ) NUadc = 0;
00889             if ( NUadc > mMaxAdc ) NUadc = mMaxAdc;
00890             if (adc > 0 && mHist[19]) mHist[19]->Fill(x->strip, NUadc);
00891             if (adc > 0 && mHist[2]) mHist[2]->Fill(NUadc);
00892         }
00893 
00897         if (mOverwrite) hit->setAdc(NUadc);
00898       } // Loop over 1 plane
00899     }
00900   }
00901   return kStOk;
00902 }
00903 
00904 //________________________________________________
00905 void StEEmcSlowMaker::setTowerGainSpread(Float_t s, Float_t mean) {
00906   LOG_INFO << "setTowerGainSpread(): gain spread: " << s << "; gain mean value: " << mean << endm;
00907   // initialize tower gain factors to mean +/- spread
00908   for ( Int_t sec=0;sec<kEEmcNumSectors;sec++ ) {
00909     for ( Int_t sub=0;sub<kEEmcNumSubSectors;sub++ ) {
00910       for ( Int_t eta=0;eta<kEEmcNumEtas;eta++ ) {
00911           Float_t f = -1.0E9;
00912           while (f <= -mean || f >= 1.0) f = gRandom->Gaus(0, s);
00913           mTowerGainFact[sec][sub][eta] = mean + f;
00914       }
00915     }
00916   }
00917 }
00918 
00919 //________________________________________________
00920 void StEEmcSlowMaker::setSmdGainSpread(Float_t s) {
00921   for (Int_t strip = 0;strip < kEEmcNumStrips;strip++) {
00922       setSmdGainSpread(s, strip);
00923   }
00924 }
00925 
00926 //________________________________________________
00927 void StEEmcSlowMaker::setSmdGainSpread(Float_t s, Int_t strip) {
00928   for (Int_t sec = 0;sec < kEEmcNumSectors;sec++) {
00929     for (Int_t uv = 0;uv < kEEmcNumSmdUVs;uv++) {
00930       setSmdGainSpread(s, sec, uv, strip);
00931     }
00932   }
00933 }
00934 
00935 //________________________________________________
00936 void StEEmcSlowMaker::setSmdGainSpread(Float_t s, Int_t sec, Int_t uv, Int_t strip) {
00937   const Float_t mean = 1.0;
00938   Float_t f = -1.0E9;
00939   while (f <= -mean || f >= 1.0) f = gRandom->Gaus(0, s);
00940   mSmdGainFact[sec][uv][strip] = mean + f;
00941 }
00942 
00943 //________________________________________________
00944 Bool_t StEEmcSlowMaker::checkDBped(const EEmcDbItem *x) {
00945   //check that channels ped and pedSigma are all >= 0
00946   return (x && (x->ped >= 0) && (x->sigPed >= 0));
00947 }
00948 
00949 //________________________________________________
00950 Float_t StEEmcSlowMaker::getPedSmear(Float_t sigPed) {
00951   //get pedestal smearing from gaussian distribution truncated at N*sigma
00952   Float_t smear = 0;
00953   if (mTruncatePedSmear > 0) {
00954     do {
00955         smear = gRandom->Gaus(0, 1.0);
00956     } while (fabs(smear) > mTruncatePedSmear);
00957   }
00958   return smear * sigPed;
00959 }
00960 
00961 //________________________________________________
00962 void StEEmcSlowMaker::setZeroAdc(StEmcCollection* emc) {
00963   StDetectorId detId[4]={kEndcapEmcTowerId,kEndcapEmcPreShowerId,kEndcapSmdUStripId,kEndcapSmdVStripId};
00964   StEmcDetector *det=0;
00965 
00966   for(int i=0; i<4; i++) {
00967     det = emc->detector(detId[i]);
00968     if (!det) continue;
00969     for (UInt_t sec = 1;sec <= det->numberOfModules();sec++) {
00970       StSPtrVecEmcRawHit &det_hits = det->module(sec)->hits();
00971       for (UInt_t ihit = 0;ihit < det_hits.size();ihit++) {
00972         StEmcRawHit *hit = det_hits[ihit];
00973         hit->setAdc(0);
00974       }
00975     }
00976   }
00977 }
00978 
00979 //________________________________________________
00980 void StEEmcSlowMaker::setZeroAdc(StMuEmcCollection *emc) {
00981   
00982   for (Int_t i = 0;i < emc->getNEndcapTowerADC();i++) 
00983     emc->setTowerADC( i+1, 0, eemc );
00984 
00985   for (Int_t i = 0;i < emc->getNEndcapPrsHits();i++) {
00986     Int_t pre,sec,eta,sub;
00987     StMuEmcHit *hit = emc->getEndcapPrsHit(i,sec,sub,eta,pre);
00988     if (hit) hit->setAdc(0);
00989   }
00990 
00991   for (Char_t uv = 'U';uv <= 'V';uv++) {
00992     Int_t sec,strip;
00993     for (Int_t i = 0;i < emc->getNEndcapSmdHits(uv);i++) {
00994       StMuEmcHit *hit=emc->getEndcapSmdHit(uv,i,sec,strip);
00995       if (hit) hit->setAdc(0);
00996     }
00997   }
00998 }
00999 
01000 //________________________________________________
01001 Float_t StEEmcSlowMaker::getMipdEdx() {
01002   // Return MIP dE/dx = 1.998 MeV/cm from the PDG book
01003   // used to simulate SMD, Pre, Post ADC response
01004   return 0.001998;
01005 }
01006 
01007 // $Log: StEEmcSlowMaker.cxx,v $
01008 // Revision 2.11  2010/09/07 22:24:52  stevens4
01009 // give access to MIP dE/dx to other makers
01010 //
01011 // Revision 2.10  2010/08/03 02:20:40  stevens4
01012 // final update from peer review
01013 //
01014 // Revision 2.9  2010/07/29 16:12:03  ogrebeny
01015 // Update after the peer review
01016 //
01017 // Revision 2.8  2010/05/03 20:47:22  ogrebeny
01018 // Some code cleanup
01019 //
01020 // Revision 2.7  2010/02/12 23:02:38  ogrebeny
01021 // By the request of the photon group, added an option to shift EEMC gains in the slow simulator.
01022 //
01023 // Revision 2.6  2009/11/23 23:44:32  ogrebeny
01024 // At Pibero's request, for the embedding infrastructure.
01025 //
01026 // Revision 2.5  2009/02/05 20:06:53  ogrebeny
01027 // Changed StEEmcDbMaker -> StEEmcDb
01028 //
01029 // Revision 2.4  2008/04/15 00:15:52  jwebb
01030 // Fixed bug in sector 12V.
01031 //
01032 // Revision 2.3  2008/04/11 14:37:17  jwebb
01033 // Added options to disable operation of individual slow simulaor subsystems.
01034 //
01035 // Revision 2.2  2007/11/28 16:17:30  jwebb
01036 // Added the following features:
01037 //
01038 // 1. User may specify the sampling fraction.
01039 //
01040 // 2. Tower and SMD gain spreads.
01041 //
01042 // Revision 2.1  2007/01/24 21:07:03  balewski
01043 // 1) no cout or printf, only new Logger
01044 // 2) EndcapMixer:
01045 //    - no assert()
01046 //    - locks out on first fatal error til the end of the job
01047 //
01048 // Revision 2.0  2007/01/13 00:03:04  jwebb
01049 // Upgrade of the slow simulator.  The following changes have been made:
01050 //
01051 // 1. Towers will always be masked out when a "fail" bit is set in the database.
01052 //    Previously this only happened if pedestals were being added, smeared.
01053 //
01054 // 2. Tower, preshower and postshower ADC values will be simulated using the
01055 //    GEANT energy loss stored in StEmcHit and StMuEmcHit.  Previously, ADC
01056 //    values from the fast simulator were used and energy loss recovered
01057 //    using gains, resulting in roundoff errors.   Note that towers still use
01058 //    the old path for MuDst-based analysis.
01059 //
01060 // 3. Tower simulation now accounts for the different light yields provided
01061 //    by the brighter scintillator and two-fiber readout in the preshower
01062 //    and postshower layers.  Previously, only the difference in thickness
01063 //    was accounted for by the GEANT simulation.
01064 //
01065 // Revision 1.5  2006/12/12 20:29:14  balewski
01066 // added hooks for Endcap embedding
01067 //
01068 // Revision 1.4  2006/08/07 18:50:11  balewski
01069 // added capabilty to run on StEvent, use se-method, see macros/ for example
01070 //
01071 // Revision 1.3  2005/09/23 17:23:01  balewski
01072 // now peds are added also to ADC of zero
01073 // fixed bug in ETOW ped smearing
01074 //
01075 // Revision 1.2  2005/09/23 01:30:10  jwebb
01076 // Tower peds now added  if option is set.
01077 //
01078 // Revision 1.1  2004/12/15 17:02:56  balewski
01079 // try 2
01080 //