CMSSW Reference Manual

00001 
00002 // File: HcalQie.cc
00003 // Description: Simulation of QIE readout for HCal
00005 
00006 #include "SimG4CMS/Calo/interface/HcalQie.h"
00007 
00008 #include "FWCore/ServiceRegistry/interface/Service.h"
00009 #include "FWCore/Utilities/interface/RandomNumberGenerator.h"
00010 #include "CLHEP/Random/RandGaussQ.h"
00011 #include "CLHEP/Random/RandPoissonQ.h"
00012 #include "FWCore/Utilities/interface/Exception.h"
00013 #include "CLHEP/Units/PhysicalConstants.h"
00014 
00015 #include <iostream>
00016 #include <iomanip>
00017 #include <cmath>
00018 
00019 HcalQie::HcalQie(edm::ParameterSet const & p) {
00020 
00021   //static SimpleConfigurable<double> p1(4.0,   "HcalQie:qToPE");
00022   //static SimpleConfigurable<int>    p2(6,     "HcalQie:BinOfMax");
00023   //static SimpleConfigurable<int>    p3(2,     "HcalQie:SignalBuckets");
00024   //static SimpleConfigurable<int>    p4(0,     "HcalQie:PreSamples");
00025   //static SimpleConfigurable<int>    p5(10,    "HcalQie:NumOfBuckets");
00026   //static SimpleConfigurable<double> p6(0.5,   "HcalQie:SigmaNoise");
00027   //static SimpleConfigurable<double> p7(0.0005,"HcalQie:EDepPerPE");
00028   //static SimpleConfigurable<int>    p8(4,     "HcalQie:BaseLine");
00029 
00030   edm::ParameterSet m_HQ  = p.getParameter<edm::ParameterSet>("HcalQie");
00031   qToPE         = m_HQ.getParameter<double>("qToPE");
00032   binOfMax      = m_HQ.getParameter<int>("BinOfMax");
00033   signalBuckets = m_HQ.getParameter<int>("SignalBuckets");
00034   preSamples    = m_HQ.getParameter<int>("PreSamples");
00035   numOfBuckets  = m_HQ.getParameter<int>("NumOfBuckets");
00036   sigma         = m_HQ.getParameter<double>("SigmaNoise");
00037   eDepPerPE     = m_HQ.getParameter<double>("EDepPerPE");
00038   int bl        = m_HQ.getParameter<int>("BaseLine");
00039 
00040   shape_  = shape();
00041   code_   = code();
00042   charge_ = charge();
00043   if (signalBuckets == 1) {
00044     phase_ = -3; rescale_ = 1.46;
00045   } else if (signalBuckets == 3) {
00046     phase_ = -1; rescale_ = 1.06;
00047   } else if (signalBuckets == 4) {
00048     phase_ =  0; rescale_ = 1.03;
00049   } else {
00050     phase_ = -2; rescale_ = 1.14; signalBuckets = 2;
00051   }
00052   weight_ = weight(binOfMax, signalBuckets, preSamples, numOfBuckets);
00053   baseline= codeToQ(bl);
00054   bmin_   = binOfMax-3;
00055   if (bmin_<0) bmin_ = 0;
00056   if (binOfMax>numOfBuckets) 
00057     bmax_ = numOfBuckets+5;
00058   else 
00059     bmax_ = binOfMax+5;
00060 
00061   edm::LogInfo("HcalSim") << "HcalQie: initialized with binOfMax " << binOfMax 
00062                           << " sample from " << bmin_ << " to " << bmax_ 
00063                           << "; signalBuckets " << signalBuckets 
00064                           << " Baseline/Phase/Scale " << baseline << "/" 
00065                           << phase_ << "/" << rescale_ << "\n"
00066                           <<"                          Noise " << sigma 
00067                           << "fC  fCToPE " << qToPE << " EDepPerPE " 
00068                           << eDepPerPE;
00069 }
00070 
00071 HcalQie::~HcalQie() {
00072   edm::LogInfo("HcalSim") << "HcalQie:: Deleting Qie";
00073 }
00074 
00075 std::vector<double> HcalQie::shape() {
00076 
00077   // pulse shape time constants in ns
00078   const float ts1  = 8.;          // scintillation time constants : 1,2,3
00079   const float ts2  = 10.;           
00080   const float ts3  = 29.3;         
00081   const float thpd = 4.;          // HPD current collection drift time
00082   const float tpre = 5.;          // preamp time constant
00083   
00084   const float wd1 = 2.;           // relative weights of decay exponents 
00085   const float wd2 = 0.7;
00086   const float wd3 = 1.;
00087 
00088   // HPD starts at I and rises to 2I in thpd of time
00089   double norm=0.0;
00090   int j, hpd_siz = (int)(thpd);
00091   std::vector<double> hpd_drift(hpd_siz);
00092   for (j=0; j<hpd_siz; j++) {
00093     double tmp = (double)j + 0.5;
00094     hpd_drift[j] = 1.0 + tmp/thpd;
00095     norm += hpd_drift[j];
00096   }
00097   // normalize integrated current to 1.0
00098   for (j=0; j<hpd_siz; j++) {
00099     hpd_drift[j] /= norm;
00100   }
00101   
00102   // Binkley shape over 6 time constants
00103   int  preamp_siz=(int)(6*tpre);
00104   std::vector<double> preamp(preamp_siz);
00105   norm = 0;
00106   for (j=0; j<preamp_siz; j++) {
00107     double tmp = (double)j + 0.5;
00108     preamp[j]  = tmp*exp(-(tmp*tmp)/(tpre*tpre));
00109     norm += preamp[j];
00110   }
00111   // normalize pulse area to 1.0
00112   for (j=0; j<preamp_siz; j++) {
00113     preamp[j] /= norm;
00114   }
00115 
00116   // ignore stochastic variation of photoelectron emission
00117   // <...>
00118   // effective tile plus wave-length shifter decay time over 4 time constants
00119 
00120   int tmax = 6 * (int)ts3;
00121   std::vector<double> scnt_decay(tmax);
00122   norm = 0;
00123   for (j=0; j<tmax; j++) {
00124     double tmp = (double)j + 0.5;
00125     scnt_decay[j] = wd1*exp(-tmp/ts1) + wd2*exp(-tmp/ts2) + wd3*exp(-tmp/ts3);
00126     norm += scnt_decay[j];
00127   }
00128   // normalize pulse area to 1.0
00129   for (j=0; j<tmax; j++) {
00130     scnt_decay[j] /= norm;
00131   }
00132 
00133   int nsiz = tmax + hpd_siz + preamp_siz + 1;
00134   std::vector<double> pulse(nsiz,0.0);  // zeroing output pulse shape
00135   norm = 0;
00136   int    i, k;
00137   for (i=0; i<tmax; i++) {
00138     int t1 = i; // and ignore jitter from optical path length
00139     for (j=0; j<hpd_siz; j++) {
00140       int t2 = t1 + j;
00141       for (k=0; k<preamp_siz; k++) {
00142         int   t3 = t2 + k;
00143         float tmp = scnt_decay[i]*hpd_drift[j]*preamp[k];
00144         pulse[t3]+= tmp;
00145         norm     += tmp;
00146       }
00147     }
00148   }
00149   
00150   // normalize for 1 GeV pulse height
00151   edm::LogInfo("HcalSim") << "HcalQie: Convoluted Shape ============== "
00152                           << "Normalisation " << norm;
00153   for (i=0; i<nsiz; i++) {
00154     pulse[i] /= norm;
00155     LogDebug("HcalSim") << "HcalQie: Pulse[" << std::setw(3) << i << "] " 
00156                         << std::setw(8) << pulse[i];
00157   }
00158   
00159   return pulse;
00160 }
00161 
00162 std::vector<int> HcalQie::code() {
00163 
00164   unsigned int CodeFADCdata[122] = {
00165       0,     1,      2,      3,      4,      5,      6,      7,      8,      9,
00166      10,    11,     12,     13,     14,     15,     16,     17,     18,     19,
00167      20,    21,     22,     23,     24,     25,     26,     27,     28,     29,
00168      30,    31,                     34,     35,     36,     37,     38,     39,
00169      40,    41,     42,     43,     44,     45,     46,     47,     48,     49,
00170      50,    51,     52,     53,     54,     55,     56,     57,     58,     59,
00171      60,    61,     62,     63,                     66,     67,     68,     69,
00172      70,    71,     72,     73,     74,     75,     76,     77,     78,     79,
00173      80,    81,     82,     83,     84,     85,     86,     87,     88,     89,
00174      90,    91,     92,     93,     94,     95,                     98,     99,
00175     100,   101,    102,    103,    104,    105,    106,    107,    108,    109,
00176     110,   111,    112,    113,    114,    115,    116,    117,    118,    119,
00177     120,   121,    122,    123,    124,    125,    126,    127 };
00178 
00179   std::vector<int> temp(122);
00180   int i;
00181   for (i = 0; i < 122; i++) 
00182     temp[i] = (int)CodeFADCdata[i];
00183 
00184   int siz = temp.size();
00185   LogDebug("HcalSim") << "HcalQie: Codes in array of size " << siz;
00186   for (i=0; i<siz; i++)
00187     LogDebug("HcalSim") << "HcalQie: Code[" << std::setw(3) << i << "] " 
00188                         << std::setw(6) << temp[i];
00189 
00190   return temp;
00191 }
00192 
00193 std::vector<double> HcalQie::charge() {
00194 
00195   double ChargeFADCdata[122] =  {
00196    -1.5,  -0.5,    0.5,    1.5,    2.5,    3.5,    4.5,    5.5,    6.5,    7.5,
00197     8.5,   9.5,   10.5,   11.5,   12.5,   13.5,   14.5,   16.5,   18.5,   20.5,
00198    22.5,  24.5,   26.5,   28.5,   31.5,   34.5,   37.5,   40.5,   44.5,   48.5,
00199    52.5,  57.5,                   62.5,   67.5,   72.5,   77.5,   82.5,   87.5,
00200    92.5,  97.5,  102.5,  107.5,  112.5,  117.5,  122.5,  127.5,  132.5,  142.5,
00201   152.5, 162.5,  172.5,  182.5,  192.5,  202.5,  217.5,  232.5,  247.5,  262.5,
00202   282.5, 302.5,  322.5,  347.5,                  372.5,  397.5,  422.5,  447.5,
00203   472.5, 497.5,  522.5,  547.5,  572.5,  597.5,  622.5,  647.5,  672.5,  697.5,
00204   722.5, 772.5,  822.5,  872.5,  922.5,  972.5, 1022.5, 1072.5, 1147.5, 1222.5,
00205  1297.5,1372.5, 1472.5, 1572.5, 1672.5, 1797.5,                 1922.5, 2047.5,
00206  2172.5,2297.5, 2422.5, 2547.5, 2672.5, 2797.5, 2922.5, 3047.5, 3172.5, 3397.5,
00207  3422.5,3547.5, 3672.5, 3922.5, 4172.5, 4422.5, 4672.5, 4922.5, 5172.5, 5422.5,
00208  5797.5,6172.5, 6547.5, 6922.5, 7422.5, 7922.5, 8422.5, 9047.5 };             
00209 
00210   std::vector<double> temp(122);
00211   int i;
00212   for (i = 0; i < 122; i++)
00213     temp[i] = (double)(ChargeFADCdata[i]);
00214 
00215   int siz = temp.size();
00216   LogDebug("HcalSim") << "HcalQie: Charges in array of size " << siz;
00217   for (i=0; i<siz; i++)
00218     LogDebug("HcalSim") << "HcalQie: Charge[" << std::setw(3) << i << "] " 
00219                         << std::setw(8) << temp[i];
00220 
00221   return temp;
00222 }
00223 
00224 std::vector<double> HcalQie::weight(int binOfMax, int mode, int npre, 
00225                                     int bucket) {
00226 
00227   std::vector<double> temp(bucket,0);
00228   int i;
00229   for (i=binOfMax-1; i<binOfMax+mode-1; i++)
00230     temp[i] = 1.;
00231   if (npre>0) {
00232     for (i=0; i<npre; i++) {
00233       int j   = binOfMax-2-i;
00234       temp[j] =-(double)mode/(double)npre;
00235     }
00236   }
00237 
00238   int siz = temp.size();
00239   LogDebug("HcalSim") << "HcalQie: Weights in array of size " << siz 
00240                       << " and Npre " << npre;
00241   for (i=0; i<siz; i++)
00242     LogDebug("HcalSim") << "HcalQie: [Weight[" << i << "] = " << temp[i];
00243 
00244   return temp;
00245 }
00246 
00247 
00248 double HcalQie::codeToQ(int ic) {
00249 
00250   double tmp=0;
00251   for (unsigned int i=0; i<code_.size(); i++) {
00252     if (ic == code_[i]) {
00253       double delta;
00254       if (i == code_.size()-1) 
00255         delta = charge_[i] - charge_[i-1];
00256       else
00257         delta = charge_[i+1] - charge_[i];
00258       tmp = charge_[i] + 0.5*delta;
00259       break;
00260     }
00261   }
00262 
00263   return tmp;
00264 }
00265 
00266 
00267 int HcalQie::getCode(double charge) {
00268 
00269   int tmp=0;
00270   for (unsigned int i=0; i<charge_.size(); i++) {
00271     if (charge < charge_[i]) {
00272       if (i>0) tmp = code_[i-1];
00273       break;
00274     }
00275   }
00276 
00277   return tmp;
00278 }
00279 
00280 
00281 double HcalQie::getShape(double time) {
00282 
00283   double tmp = 0;
00284   int    k = (int) (time + 0.5);
00285   if (k>=0 && k<((int)(shape_.size())-1))
00286     tmp = 0.5*(shape_[k]+shape_[k+1]);
00287   
00288   return tmp;
00289 }
00290 
00291 
00292 std::vector<int> HcalQie::getCode(int nht, std::vector<CaloHit> hitbuf) {
00293 
00294   const double  bunchSpace=25.;
00295   int nmax = (bmax_ > numOfBuckets ? bmax_ : numOfBuckets);
00296   std::vector<double> work(nmax);
00297 
00298   edm::Service<edm::RandomNumberGenerator> rng;
00299   if ( ! rng.isAvailable()) {
00300     throw cms::Exception("Configuration")
00301       << "HcalQIE requires the RandomNumberGeneratorService\n"
00302       << "which is not present in the configuration file. "
00303       << "You must add the service\n in the configuration file or "
00304       << "remove the modules that require it.";
00305   }
00306   CLHEP::RandGaussQ  randGauss(rng->getEngine(), baseline,sigma);
00307   CLHEP::RandPoissonQ randPoisson(rng->getEngine());
00308 
00309 
00310   // Noise in the channel
00311   for (int i=0; i<numOfBuckets; i++) 
00312     work[i] = randGauss.fire();
00313 
00314   LogDebug("HcalSim") << "HcalQie::getCode: Noise with baseline " << baseline 
00315                       << " width " << sigma << " and " << nht << " hits";
00316   for (int i=0; i<numOfBuckets; i++) 
00317     LogDebug("HcalSim") << "HcalQie: Code[" << i << "] = " << work[i];
00318   
00319   double etot=0, esum=0, photons=0;
00320   if (nht>0) {
00321   
00322     // Sort the hits
00323     std::vector<CaloHit*> hits(nht);
00324     std::vector<CaloHit*>::iterator k1, k2;
00325     int kk;
00326     for (kk = 0; kk < nht; kk++) {
00327       hits[kk] = &hitbuf[kk];
00328     }
00329     sort(hits.begin(),hits.end(),CaloHitMore());
00330     
00331     // Energy deposits
00332     for (kk = 0, k1 = hits.begin(); k1 != hits.end(); kk++, k1++) {
00333       double ehit  = (**k1).e();
00334       double jitter= (**k1).t();
00335       int    jump  = 0;
00336       for (k2 = k1+1; k2 != hits.end() && (jitter-(**k2).t())<1. &&
00337              (jitter-(**k2).t())>-1.; k2++) {
00338         ehit += (**k2).e();
00339         jump++;
00340       }
00341 
00342       double avpe  = ehit/eDepPerPE;
00343       double photo = randPoisson.fire(avpe);
00344       etot   += ehit;
00345       photons+= photo; 
00346       LogDebug("HcalSim") << "HcalQie::getCode: Hit " << kk << ":" << kk+jump 
00347                           << " Energy deposit " << ehit << " Time " << jitter
00348                           << " Average and true no of PE " << avpe << " " 
00349                           << photo;
00350 
00351       double bintime = jitter - phase_ - bunchSpace*(binOfMax-bmin_);
00352       LogDebug("HcalSim") << "HcalQie::getCode: phase " << phase_  
00353                           << " binTime " << bintime;
00354       std::vector<double> binsum(nmax,0);
00355       double norm=0, sum=0.;
00356       for (int i=bmin_; i<bmax_; i++) {
00357         bintime += bunchSpace;
00358         for (int j=0; j<(int)(bunchSpace); j++) {
00359           double tim = bintime + j;
00360           double tmp = getShape(tim);
00361           binsum[i] += tmp;
00362         }
00363         sum += binsum[i];
00364       }
00365 
00366       if (sum>0) norm = (photo/(sum*qToPE));
00367       LogDebug("HcalSim") << "HcalQie::getCode: PE " << photo << " Sum " << sum
00368                           << " Norm. " << norm;
00369       for (int i=bmin_; i<bmax_; i++)
00370         work[i] += binsum[i]*norm;
00371 
00372       kk += jump;
00373       k1 += jump;
00374     }
00375   }
00376 
00377   std::vector<int> temp(numOfBuckets,0);
00378   for (int i=0; i<numOfBuckets; i++) {
00379     temp[i] = getCode(work[i]);
00380     esum   += work[i];
00381   }
00382   LogDebug("HcalSim") << "HcalQie::getCode: Input " << etot << " GeV; Photons "
00383                       << photons << ";  Output " << esum << " fc";
00384 
00385   return temp;
00386 }
00387 
00388 
00389 double HcalQie::getEnergy(std::vector<int> code) {
00390 
00391   std::vector<double> work(numOfBuckets);
00392   double sum=0;
00393   for (int i=0; i<numOfBuckets; i++) {
00394     work[i] = codeToQ (code[i]);
00395     sum += work[i]*weight_[i];
00396     LogDebug("HcalSim") << "HcalQie::getEnergy: " << i << " code " << code[i] 
00397                         << " PE " << work[i];
00398   }
00399 
00400   double tmp;
00401   if (preSamples == 0) {
00402     tmp = (sum - signalBuckets*baseline) * rescale_ * qToPE * eDepPerPE;
00403   } else {
00404     tmp = sum * rescale_ * qToPE;
00405   }
00406   LogDebug("HcalSim") << "HcalQie::getEnergy: PE " << sum*qToPE << " Energy " 
00407                       << tmp << " GeV";
00408 
00409   return tmp;
00410 }