Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 12:29:53 

0001 ///////////////////////////////////////////////////////////////////////////////
0002 // File: HcalQie.cc
0003 // Description: Simulation of QIE readout for HCal
0004 ///////////////////////////////////////////////////////////////////////////////
0005 
0006 #include "SimG4CMS/Calo/interface/HcalQie.h"
0007 
0008 #include "CLHEP/Random/RandGaussQ.h"
0009 #include "CLHEP/Random/RandPoissonQ.h"
0010 #include "CLHEP/Units/GlobalPhysicalConstants.h"
0011 
0012 #include <iostream>
0013 #include <iomanip>
0014 #include <cmath>
0015 
0016 //#define EDM_ML_DEBUG
0017 
0018 HcalQie::HcalQie(edm::ParameterSet const& p) {
0019   //static SimpleConfigurable<double> p1(4.0,   "HcalQie:qToPE");
0020   //static SimpleConfigurable<int>    p2(6,     "HcalQie:BinOfMax");
0021   //static SimpleConfigurable<int>    p3(2,     "HcalQie:SignalBuckets");
0022   //static SimpleConfigurable<int>    p4(0,     "HcalQie:PreSamples");
0023   //static SimpleConfigurable<int>    p5(10,    "HcalQie:NumOfBuckets");
0024   //static SimpleConfigurable<double> p6(0.5,   "HcalQie:SigmaNoise");
0025   //static SimpleConfigurable<double> p7(0.0005,"HcalQie:EDepPerPE");
0026   //static SimpleConfigurable<int>    p8(4,     "HcalQie:BaseLine");
0027 
0028   edm::ParameterSet m_HQ = p.getParameter<edm::ParameterSet>("HcalQie");
0029   qToPE = m_HQ.getParameter<double>("qToPE");
0030   binOfMax = m_HQ.getParameter<int>("BinOfMax");
0031   signalBuckets = m_HQ.getParameter<int>("SignalBuckets");
0032   preSamples = m_HQ.getParameter<int>("PreSamples");
0033   numOfBuckets = m_HQ.getParameter<int>("NumOfBuckets");
0034   sigma = m_HQ.getParameter<double>("SigmaNoise");
0035   eDepPerPE = m_HQ.getParameter<double>("EDepPerPE");
0036   int bl = m_HQ.getParameter<int>("BaseLine");
0037 
0038   shape_ = shape();
0039   code_ = code();
0040   charge_ = charge();
0041   if (signalBuckets == 1) {
0042     phase_ = -3;
0043     rescale_ = 1.46;
0044   } else if (signalBuckets == 3) {
0045     phase_ = -1;
0046     rescale_ = 1.06;
0047   } else if (signalBuckets == 4) {
0048     phase_ = 0;
0049     rescale_ = 1.03;
0050   } else {
0051     phase_ = -2;
0052     rescale_ = 1.14;
0053     signalBuckets = 2;
0054   }
0055   weight_ = weight(binOfMax, signalBuckets, preSamples, numOfBuckets);
0056   baseline = codeToQ(bl);
0057   bmin_ = binOfMax - 3;
0058   if (bmin_ < 0)
0059     bmin_ = 0;
0060   if (binOfMax > numOfBuckets)
0061     bmax_ = numOfBuckets + 5;
0062   else
0063     bmax_ = binOfMax + 5;
0064 
0065   edm::LogVerbatim("HcalSim") << "HcalQie: initialized with binOfMax " << binOfMax << " sample from " << bmin_ << " to "
0066                               << bmax_ << "; signalBuckets " << signalBuckets << " Baseline/Phase/Scale " << baseline
0067                               << "/" << phase_ << "/" << rescale_ << "\n                          Noise " << sigma
0068                               << "fC  fCToPE " << qToPE << " EDepPerPE " << eDepPerPE;
0069 }
0070 
0071 HcalQie::~HcalQie() { edm::LogVerbatim("HcalSim") << "HcalQie:: Deleting Qie"; }
0072 
0073 std::vector<double> HcalQie::shape() {
0074   // pulse shape time constants in ns
0075   const float ts1 = 8.;  // scintillation time constants : 1,2,3
0076   const float ts2 = 10.;
0077   const float ts3 = 29.3;
0078   const float thpd = 4.;  // HPD current collection drift time
0079   const float tpre = 5.;  // preamp time constant
0080 
0081   const float wd1 = 2.;  // relative weights of decay exponents
0082   const float wd2 = 0.7;
0083   const float wd3 = 1.;
0084 
0085   // HPD starts at I and rises to 2I in thpd of time
0086   double norm = 0.0;
0087   int j, hpd_siz = (int)(thpd);
0088   std::vector<double> hpd_drift(hpd_siz);
0089   for (j = 0; j < hpd_siz; j++) {
0090     double tmp = (double)j + 0.5;
0091     hpd_drift[j] = 1.0 + tmp / thpd;
0092     norm += hpd_drift[j];
0093   }
0094   // normalize integrated current to 1.0
0095   for (j = 0; j < hpd_siz; j++) {
0096     hpd_drift[j] /= norm;
0097   }
0098 
0099   // Binkley shape over 6 time constants
0100   int preamp_siz = (int)(6 * tpre);
0101   std::vector<double> preamp(preamp_siz);
0102   norm = 0;
0103   for (j = 0; j < preamp_siz; j++) {
0104     double tmp = (double)j + 0.5;
0105     preamp[j] = tmp * exp(-(tmp * tmp) / (tpre * tpre));
0106     norm += preamp[j];
0107   }
0108   // normalize pulse area to 1.0
0109   for (j = 0; j < preamp_siz; j++) {
0110     preamp[j] /= norm;
0111   }
0112 
0113   // ignore stochastic variation of photoelectron emission
0114   // <...>
0115   // effective tile plus wave-length shifter decay time over 4 time constants
0116 
0117   int tmax = 6 * (int)ts3;
0118   std::vector<double> scnt_decay(tmax);
0119   norm = 0;
0120   for (j = 0; j < tmax; j++) {
0121     double tmp = (double)j + 0.5;
0122     scnt_decay[j] = wd1 * exp(-tmp / ts1) + wd2 * exp(-tmp / ts2) + wd3 * exp(-tmp / ts3);
0123     norm += scnt_decay[j];
0124   }
0125   // normalize pulse area to 1.0
0126   for (j = 0; j < tmax; j++) {
0127     scnt_decay[j] /= norm;
0128   }
0129 
0130   int nsiz = tmax + hpd_siz + preamp_siz + 1;
0131   std::vector<double> pulse(nsiz, 0.0);  // zeroing output pulse shape
0132   norm = 0;
0133   int i, k;
0134   for (i = 0; i < tmax; i++) {
0135     int t1 = i;  // and ignore jitter from optical path length
0136     for (j = 0; j < hpd_siz; j++) {
0137       int t2 = t1 + j;
0138       for (k = 0; k < preamp_siz; k++) {
0139         int t3 = t2 + k;
0140         float tmp = scnt_decay[i] * hpd_drift[j] * preamp[k];
0141         pulse[t3] += tmp;
0142         norm += tmp;
0143       }
0144     }
0145   }
0146 
0147   // normalize for 1 GeV pulse height
0148   edm::LogVerbatim("HcalSim") << "HcalQie: Convoluted Shape ============== Normalisation " << norm;
0149   for (i = 0; i < nsiz; i++) {
0150     pulse[i] /= norm;
0151 #ifdef EDM_ML_DEBUG
0152     edm::LogVerbatim("HcalSim") << "HcalQie: Pulse[" << std::setw(3) << i << "] " << std::setw(8) << pulse[i];
0153 #endif
0154   }
0155 
0156   return pulse;
0157 }
0158 
0159 std::vector<int> HcalQie::code() {
0160   unsigned int CodeFADCdata[122] = {
0161       0,   1,   2,   3,   4,   5,   6,   7,   8,   9,   10,  11,  12,  13,  14,  15,  16,  17,  18,  19,  20,
0162       21,  22,  23,  24,  25,  26,  27,  28,  29,  30,  31,  34,  35,  36,  37,  38,  39,  40,  41,  42,  43,
0163       44,  45,  46,  47,  48,  49,  50,  51,  52,  53,  54,  55,  56,  57,  58,  59,  60,  61,  62,  63,  66,
0164       67,  68,  69,  70,  71,  72,  73,  74,  75,  76,  77,  78,  79,  80,  81,  82,  83,  84,  85,  86,  87,
0165       88,  89,  90,  91,  92,  93,  94,  95,  98,  99,  100, 101, 102, 103, 104, 105, 106, 107, 108, 109, 110,
0166       111, 112, 113, 114, 115, 116, 117, 118, 119, 120, 121, 122, 123, 124, 125, 126, 127};
0167 
0168   std::vector<int> temp(122);
0169   int i;
0170   for (i = 0; i < 122; i++)
0171     temp[i] = (int)CodeFADCdata[i];
0172 
0173 #ifdef EDM_ML_DEBUG
0174   int siz = temp.size();
0175   edm::LogVerbatim("HcalSim") << "HcalQie: Codes in array of size " << siz;
0176   for (i = 0; i < siz; i++)
0177     edm::LogVerbatim("HcalSim") << "HcalQie: Code[" << std::setw(3) << i << "] " << std::setw(6) << temp[i];
0178 #endif
0179   return temp;
0180 }
0181 
0182 std::vector<double> HcalQie::charge() {
0183   double ChargeFADCdata[122] = {
0184       -1.5,   -0.5,   0.5,    1.5,    2.5,    3.5,    4.5,    5.5,    6.5,    7.5,    8.5,    9.5,    10.5,   11.5,
0185       12.5,   13.5,   14.5,   16.5,   18.5,   20.5,   22.5,   24.5,   26.5,   28.5,   31.5,   34.5,   37.5,   40.5,
0186       44.5,   48.5,   52.5,   57.5,   62.5,   67.5,   72.5,   77.5,   82.5,   87.5,   92.5,   97.5,   102.5,  107.5,
0187       112.5,  117.5,  122.5,  127.5,  132.5,  142.5,  152.5,  162.5,  172.5,  182.5,  192.5,  202.5,  217.5,  232.5,
0188       247.5,  262.5,  282.5,  302.5,  322.5,  347.5,  372.5,  397.5,  422.5,  447.5,  472.5,  497.5,  522.5,  547.5,
0189       572.5,  597.5,  622.5,  647.5,  672.5,  697.5,  722.5,  772.5,  822.5,  872.5,  922.5,  972.5,  1022.5, 1072.5,
0190       1147.5, 1222.5, 1297.5, 1372.5, 1472.5, 1572.5, 1672.5, 1797.5, 1922.5, 2047.5, 2172.5, 2297.5, 2422.5, 2547.5,
0191       2672.5, 2797.5, 2922.5, 3047.5, 3172.5, 3397.5, 3422.5, 3547.5, 3672.5, 3922.5, 4172.5, 4422.5, 4672.5, 4922.5,
0192       5172.5, 5422.5, 5797.5, 6172.5, 6547.5, 6922.5, 7422.5, 7922.5, 8422.5, 9047.5};
0193 
0194   std::vector<double> temp(122);
0195   int i;
0196   for (i = 0; i < 122; i++)
0197     temp[i] = (double)(ChargeFADCdata[i]);
0198 
0199 #ifdef EDM_ML_DEBUG
0200   int siz = temp.size();
0201   edm::LogVerbatim("HcalSim") << "HcalQie: Charges in array of size " << siz;
0202   for (i = 0; i < siz; i++)
0203     edm::LogVerbatim("HcalSim") << "HcalQie: Charge[" << std::setw(3) << i << "] " << std::setw(8) << temp[i];
0204 #endif
0205   return temp;
0206 }
0207 
0208 std::vector<double> HcalQie::weight(int binOfMax, int mode, int npre, int bucket) {
0209   std::vector<double> temp(bucket, 0);
0210   int i;
0211   for (i = binOfMax - 1; i < binOfMax + mode - 1; i++)
0212     temp[i] = 1.;
0213   if (npre > 0) {
0214     for (i = 0; i < npre; i++) {
0215       int j = binOfMax - 2 - i;
0216       temp[j] = -(double)mode / (double)npre;
0217     }
0218   }
0219 
0220 #ifdef EDM_ML_DEBUG
0221   int siz = temp.size();
0222   edm::LogVerbatim("HcalSim") << "HcalQie: Weights in array of size " << siz << " and Npre " << npre;
0223   for (i = 0; i < siz; i++)
0224     edm::LogVerbatim("HcalSim") << "HcalQie: [Weight[" << i << "] = " << temp[i];
0225 #endif
0226   return temp;
0227 }
0228 
0229 double HcalQie::codeToQ(int ic) {
0230   double tmp = 0;
0231   for (unsigned int i = 0; i < code_.size(); i++) {
0232     if (ic == code_[i]) {
0233       double delta;
0234       if (i == code_.size() - 1)
0235         delta = charge_[i] - charge_[i - 1];
0236       else
0237         delta = charge_[i + 1] - charge_[i];
0238       tmp = charge_[i] + 0.5 * delta;
0239       break;
0240     }
0241   }
0242 
0243   return tmp;
0244 }
0245 
0246 int HcalQie::getCode(double charge) {
0247   int tmp = 0;
0248   for (unsigned int i = 0; i < charge_.size(); i++) {
0249     if (charge < charge_[i]) {
0250       if (i > 0)
0251         tmp = code_[i - 1];
0252       break;
0253     }
0254   }
0255 
0256   return tmp;
0257 }
0258 
0259 double HcalQie::getShape(double time) {
0260   double tmp = 0;
0261   int k = (int)(time + 0.5);
0262   if (k >= 0 && k < ((int)(shape_.size()) - 1))
0263     tmp = 0.5 * (shape_[k] + shape_[k + 1]);
0264 
0265   return tmp;
0266 }
0267 
0268 std::vector<int> HcalQie::getCode(int nht, const std::vector<CaloHit>& hitbuf, CLHEP::HepRandomEngine* engine) {
0269   const double bunchSpace = 25.;
0270   int nmax = (bmax_ > numOfBuckets ? bmax_ : numOfBuckets);
0271   std::vector<double> work(nmax);
0272 
0273   // Noise in the channel
0274   for (int i = 0; i < numOfBuckets; i++)
0275     work[i] = CLHEP::RandGaussQ::shoot(engine, baseline, sigma);
0276 
0277 #ifdef EDM_ML_DEBUG
0278   edm::LogVerbatim("HcalSim") << "HcalQie::getCode: Noise with baseline " << baseline << " width " << sigma << " and "
0279                               << nht << " hits";
0280   for (int i = 0; i < numOfBuckets; i++)
0281     edm::LogVerbatim("HcalSim") << "HcalQie: Code[" << i << "] = " << work[i];
0282   double etot = 0, esum = 0, photons = 0;
0283 #endif
0284   if (nht > 0) {
0285     // Sort the hits
0286     std::vector<const CaloHit*> hits(nht);
0287     std::vector<const CaloHit*>::iterator k1, k2;
0288     int kk;
0289     for (kk = 0; kk < nht; kk++) {
0290       hits[kk] = &hitbuf[kk];
0291     }
0292     sort(hits.begin(), hits.end(), CaloHitMore());
0293 
0294     // Energy deposits
0295     for (kk = 0, k1 = hits.begin(); k1 != hits.end(); kk++, k1++) {
0296       double ehit = (**k1).e();
0297       double jitter = (**k1).t();
0298       int jump = 0;
0299       for (k2 = k1 + 1; k2 != hits.end() && (jitter - (**k2).t()) < 1. && (jitter - (**k2).t()) > -1.; k2++) {
0300         ehit += (**k2).e();
0301         jump++;
0302       }
0303 
0304       double avpe = ehit / eDepPerPE;
0305       CLHEP::RandPoissonQ randPoissonQ(*engine, avpe);
0306       double photo = randPoissonQ.fire();
0307 #ifdef EDM_ML_DEBUG
0308       etot += ehit;
0309       photons += photo;
0310       edm::LogVerbatim("HcalSim") << "HcalQie::getCode: Hit " << kk << ":" << kk + jump << " Energy deposit " << ehit
0311                                   << " Time " << jitter << " Average and true no of PE " << avpe << " " << photo;
0312 #endif
0313       double bintime = jitter - phase_ - bunchSpace * (binOfMax - bmin_);
0314 #ifdef EDM_ML_DEBUG
0315       edm::LogVerbatim("HcalSim") << "HcalQie::getCode: phase " << phase_ << " binTime " << bintime;
0316 #endif
0317       std::vector<double> binsum(nmax, 0);
0318       double norm = 0, sum = 0.;
0319       for (int i = bmin_; i < bmax_; i++) {
0320         bintime += bunchSpace;
0321         for (int j = 0; j < (int)(bunchSpace); j++) {
0322           double tim = bintime + j;
0323           double tmp = getShape(tim);
0324           binsum[i] += tmp;
0325         }
0326         sum += binsum[i];
0327       }
0328 
0329       if (sum > 0)
0330         norm = (photo / (sum * qToPE));
0331 #ifdef EDM_ML_DEBUG
0332       edm::LogVerbatim("HcalSim") << "HcalQie::getCode: PE " << photo << " Sum " << sum << " Norm. " << norm;
0333 #endif
0334       for (int i = bmin_; i < bmax_; i++)
0335         work[i] += binsum[i] * norm;
0336 
0337       kk += jump;
0338       k1 += jump;
0339     }
0340   }
0341 
0342   std::vector<int> temp(numOfBuckets, 0);
0343   for (int i = 0; i < numOfBuckets; i++) {
0344     temp[i] = getCode(work[i]);
0345 #ifdef EDM_ML_DEBUG
0346     esum += work[i];
0347 #endif
0348   }
0349 #ifdef EDM_ML_DEBUG
0350   edm::LogVerbatim("HcalSim") << "HcalQie::getCode: Input " << etot << " GeV; Photons " << photons << ";  Output "
0351                               << esum << " fc";
0352 #endif
0353   return temp;
0354 }
0355 
0356 double HcalQie::getEnergy(const std::vector<int>& code) {
0357   std::vector<double> work(numOfBuckets);
0358   double sum = 0;
0359   for (int i = 0; i < numOfBuckets; i++) {
0360     work[i] = codeToQ(code[i]);
0361     sum += work[i] * weight_[i];
0362 #ifdef EDM_ML_DEBUG
0363     edm::LogVerbatim("HcalSim") << "HcalQie::getEnergy: " << i << " code " << code[i] << " PE " << work[i];
0364 #endif
0365   }
0366 
0367   double tmp;
0368   if (preSamples == 0) {
0369     tmp = (sum - signalBuckets * baseline) * rescale_ * qToPE * eDepPerPE;
0370   } else {
0371     tmp = sum * rescale_ * qToPE;
0372   }
0373 #ifdef EDM_ML_DEBUG
0374   edm::LogVerbatim("HcalSim") << "HcalQie::getEnergy: PE " << sum * qToPE << " Energy " << tmp << " GeV";
0375 #endif
0376   return tmp;
0377 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!