#include <HcalQie.h>

Public Member Functions
std::vector< int >	getCode (int, std::vector< CaloHit >)
double	getEnergy (std::vector< int >)
	HcalQie (edm::ParameterSet const &p)
virtual	~HcalQie ()
Private Member Functions
std::vector< double >	charge ()
std::vector< int >	code ()
double	codeToQ (int ic)
int	getCode (double charge)
double	getShape (double time)
std::vector< double >	shape ()
std::vector< double >	weight (int binofmax, int mode, int npre, int numbucket)
Private Attributes
double	baseline
int	binOfMax
int	bmax_
int	bmin_
std::vector< double >	charge_
std::vector< int >	code_
double	eDepPerPE
int	numOfBuckets
double	phase_
int	preSamples
double	qToPE
double	rescale_
std::vector< double >	shape_
double	sigma
int	signalBuckets
std::vector< double >	weight_

Detailed Description

Definition at line 14 of file HcalQie.h.

Constructor & Destructor Documentation

HcalQie::HcalQie ( edm::ParameterSet const & p )

Definition at line 21 of file HcalQie.cc.

References baseline, binOfMax, bmax_, bmin_, charge(), charge_, code(), code_, codeToQ(), eDepPerPE, edm::ParameterSet::getParameter(), numOfBuckets, phase_, preSamples, qToPE, rescale_, shape(), shape_, sigma, signalBuckets, weight(), and weight_.

                                          {

  //static SimpleConfigurable<double> p1(4.0,   "HcalQie:qToPE");
  //static SimpleConfigurable<int>    p2(6,     "HcalQie:BinOfMax");
  //static SimpleConfigurable<int>    p3(2,     "HcalQie:SignalBuckets");
  //static SimpleConfigurable<int>    p4(0,     "HcalQie:PreSamples");
  //static SimpleConfigurable<int>    p5(10,    "HcalQie:NumOfBuckets");
  //static SimpleConfigurable<double> p6(0.5,   "HcalQie:SigmaNoise");
  //static SimpleConfigurable<double> p7(0.0005,"HcalQie:EDepPerPE");
  //static SimpleConfigurable<int>    p8(4,     "HcalQie:BaseLine");

  edm::ParameterSet m_HQ  = p.getParameter<edm::ParameterSet>("HcalQie");
  qToPE         = m_HQ.getParameter<double>("qToPE");
  binOfMax      = m_HQ.getParameter<int>("BinOfMax");
  signalBuckets = m_HQ.getParameter<int>("SignalBuckets");
  preSamples    = m_HQ.getParameter<int>("PreSamples");
  numOfBuckets  = m_HQ.getParameter<int>("NumOfBuckets");
  sigma         = m_HQ.getParameter<double>("SigmaNoise");
  eDepPerPE     = m_HQ.getParameter<double>("EDepPerPE");
  int bl        = m_HQ.getParameter<int>("BaseLine");

  shape_  = shape();
  code_   = code();
  charge_ = charge();
  if (signalBuckets == 1) {
    phase_ = -3; rescale_ = 1.46;
  } else if (signalBuckets == 3) {
    phase_ = -1; rescale_ = 1.06;
  } else if (signalBuckets == 4) {
    phase_ =  0; rescale_ = 1.03;
  } else {
    phase_ = -2; rescale_ = 1.14; signalBuckets = 2;
  }
  weight_ = weight(binOfMax, signalBuckets, preSamples, numOfBuckets);
  baseline= codeToQ(bl);
  bmin_   = binOfMax-3;
  if (bmin_<0) bmin_ = 0;
  if (binOfMax>numOfBuckets) 
    bmax_ = numOfBuckets+5;
  else 
    bmax_ = binOfMax+5;

  edm::LogInfo("HcalSim") << "HcalQie: initialized with binOfMax " << binOfMax 
                          << " sample from " << bmin_ << " to " << bmax_ 
                          << "; signalBuckets " << signalBuckets 
                          << " Baseline/Phase/Scale " << baseline << "/" 
                          << phase_ << "/" << rescale_ << "\n"
                          <<"                          Noise " << sigma 
                          << "fC  fCToPE " << qToPE << " EDepPerPE " 
                          << eDepPerPE;
}

HcalQie::~HcalQie ( ) [virtual]

Definition at line 73 of file HcalQie.cc.

                  {
  edm::LogInfo("HcalSim") << "HcalQie:: Deleting Qie";
}

Member Function Documentation

std::vector< double > HcalQie::charge ( void ) [private]

Definition at line 198 of file HcalQie.cc.

References i, LogDebug, and groupFilesInBlocks::temp.

Referenced by HcalQie().

                                  {

  double ChargeFADCdata[122] =  {
   -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,   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,   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,  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,  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,  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, 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, 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, 5172.5, 5422.5,
 5797.5,6172.5, 6547.5, 6922.5, 7422.5, 7922.5, 8422.5, 9047.5 };             

  std::vector<double> temp(122);
  int i;
  for (i = 0; i < 122; i++)
    temp[i] = (double)(ChargeFADCdata[i]);

#ifdef DebugLog
  int siz = temp.size();
  LogDebug("HcalSim") << "HcalQie: Charges in array of size " << siz;
  for (i=0; i<siz; i++)
    LogDebug("HcalSim") << "HcalQie: Charge[" << std::setw(3) << i << "] " 
                        << std::setw(8) << temp[i];
#endif
  return temp;
}

std::vector< int > HcalQie::code ( ) [private]

Definition at line 166 of file HcalQie.cc.

References i, LogDebug, and groupFilesInBlocks::temp.

Referenced by HcalQie().

                             {

  unsigned int CodeFADCdata[122] = {
      0,     1,      2,      3,      4,      5,      6,      7,      8,      9,
     10,    11,     12,     13,     14,     15,     16,     17,     18,     19,
     20,    21,     22,     23,     24,     25,     26,     27,     28,     29,
     30,    31,                     34,     35,     36,     37,     38,     39,
     40,    41,     42,     43,     44,     45,     46,     47,     48,     49,
     50,    51,     52,     53,     54,     55,     56,     57,     58,     59,
     60,    61,     62,     63,                     66,     67,     68,     69,
     70,    71,     72,     73,     74,     75,     76,     77,     78,     79,
     80,    81,     82,     83,     84,     85,     86,     87,     88,     89,
     90,    91,     92,     93,     94,     95,                     98,     99,
    100,   101,    102,    103,    104,    105,    106,    107,    108,    109,
    110,   111,    112,    113,    114,    115,    116,    117,    118,    119,
    120,   121,    122,    123,    124,    125,    126,    127 };

  std::vector<int> temp(122);
  int i;
  for (i = 0; i < 122; i++) 
    temp[i] = (int)CodeFADCdata[i];

#ifdef DebugLog
  int siz = temp.size();
  LogDebug("HcalSim") << "HcalQie: Codes in array of size " << siz;
  for (i=0; i<siz; i++)
    LogDebug("HcalSim") << "HcalQie: Code[" << std::setw(3) << i << "] " 
                        << std::setw(6) << temp[i];
#endif
  return temp;
}

double HcalQie::codeToQ ( int ic ) [private]

Definition at line 255 of file HcalQie.cc.

References charge_, code_, delta, i, and tmp.

Referenced by getEnergy(), and HcalQie().

                              {

  double tmp=0;
  for (unsigned int i=0; i<code_.size(); i++) {
    if (ic == code_[i]) {
      double delta;
      if (i == code_.size()-1) 
        delta = charge_[i] - charge_[i-1];
      else
        delta = charge_[i+1] - charge_[i];
      tmp = charge_[i] + 0.5*delta;
      break;
    }
  }

  return tmp;
}

int HcalQie::getCode ( double charge ) [private]

Definition at line 274 of file HcalQie.cc.

References charge_, code_, i, and tmp.

                                  {

  int tmp=0;
  for (unsigned int i=0; i<charge_.size(); i++) {
    if (charge < charge_[i]) {
      if (i>0) tmp = code_[i-1];
      break;
    }
  }

  return tmp;
}

std::vector< int > HcalQie::getCode	(	int	nht,
		std::vector< CaloHit >	hitbuf
	)

Definition at line 299 of file HcalQie.cc.

References baseline, binOfMax, bmax_, bmin_, eDepPerPE, Exception, getShape(), i, edm::Service< T >::isAvailable(), j, LogDebug, numOfBuckets, phase_, interactiveExample::photons, qToPE, sigma, python::multivaluedict::sort(), groupFilesInBlocks::temp, and tmp.

Referenced by HcalTB04Analysis::qieAnalysis(), and HcalTestAnalysis::qieAnalysis().

                                                                  {

  const double  bunchSpace=25.;
  int nmax = (bmax_ > numOfBuckets ? bmax_ : numOfBuckets);
  std::vector<double> work(nmax);

  edm::Service<edm::RandomNumberGenerator> rng;
  if ( ! rng.isAvailable()) {
    throw cms::Exception("Configuration")
      << "HcalQIE requires the RandomNumberGeneratorService\n"
      << "which is not present in the configuration file. "
      << "You must add the service\n in the configuration file or "
      << "remove the modules that require it.";
  }
  CLHEP::RandGaussQ  randGauss(rng->getEngine(), baseline,sigma);
  CLHEP::RandPoissonQ randPoisson(rng->getEngine());


  // Noise in the channel
  for (int i=0; i<numOfBuckets; i++) 
    work[i] = randGauss.fire();

#ifdef DebugLog
  LogDebug("HcalSim") << "HcalQie::getCode: Noise with baseline " << baseline 
                      << " width " << sigma << " and " << nht << " hits";
  for (int i=0; i<numOfBuckets; i++) 
    LogDebug("HcalSim") << "HcalQie: Code[" << i << "] = " << work[i];
#endif  
  double etot=0, esum=0, photons=0;
  if (nht>0) {
  
    // Sort the hits
    std::vector<CaloHit*> hits(nht);
    std::vector<CaloHit*>::iterator k1, k2;
    int kk;
    for (kk = 0; kk < nht; kk++) {
      hits[kk] = &hitbuf[kk];
    }
    sort(hits.begin(),hits.end(),CaloHitMore());
    
    // Energy deposits
    for (kk = 0, k1 = hits.begin(); k1 != hits.end(); kk++, k1++) {
      double ehit  = (**k1).e();
      double jitter= (**k1).t();
      int    jump  = 0;
      for (k2 = k1+1; k2 != hits.end() && (jitter-(**k2).t())<1. &&
             (jitter-(**k2).t())>-1.; k2++) {
        ehit += (**k2).e();
        jump++;
      }

      double avpe  = ehit/eDepPerPE;
      double photo = randPoisson.fire(avpe);
      etot   += ehit;
      photons+= photo; 
#ifdef DebugLog
      LogDebug("HcalSim") << "HcalQie::getCode: Hit " << kk << ":" << kk+jump 
                          << " Energy deposit " << ehit << " Time " << jitter
                          << " Average and true no of PE " << avpe << " " 
                          << photo;
#endif
      double bintime = jitter - phase_ - bunchSpace*(binOfMax-bmin_);
#ifdef DebugLog
      LogDebug("HcalSim") << "HcalQie::getCode: phase " << phase_  
                          << " binTime " << bintime;
#endif
      std::vector<double> binsum(nmax,0);
      double norm=0, sum=0.;
      for (int i=bmin_; i<bmax_; i++) {
        bintime += bunchSpace;
        for (int j=0; j<(int)(bunchSpace); j++) {
          double tim = bintime + j;
          double tmp = getShape(tim);
          binsum[i] += tmp;
        }
        sum += binsum[i];
      }

      if (sum>0) norm = (photo/(sum*qToPE));
#ifdef DebugLog
      LogDebug("HcalSim") << "HcalQie::getCode: PE " << photo << " Sum " << sum
                          << " Norm. " << norm;
#endif
      for (int i=bmin_; i<bmax_; i++)
        work[i] += binsum[i]*norm;

      kk += jump;
      k1 += jump;
    }
  }

  std::vector<int> temp(numOfBuckets,0);
  for (int i=0; i<numOfBuckets; i++) {
    temp[i] = getCode(work[i]);
    esum   += work[i];
  }
#ifdef DebugLog
  LogDebug("HcalSim") << "HcalQie::getCode: Input " << etot << " GeV; Photons "
                      << photons << ";  Output " << esum << " fc";
#endif
  return temp;
}

double HcalQie::getEnergy ( std::vector< int > code )

Definition at line 403 of file HcalQie.cc.

References baseline, codeToQ(), eDepPerPE, i, LogDebug, numOfBuckets, preSamples, qToPE, rescale_, signalBuckets, tmp, and weight_.

Referenced by HcalTB04Analysis::qieAnalysis(), and HcalTestAnalysis::qieAnalysis().

                                             {

  std::vector<double> work(numOfBuckets);
  double sum=0;
  for (int i=0; i<numOfBuckets; i++) {
    work[i] = codeToQ (code[i]);
    sum += work[i]*weight_[i];
#ifdef DebugLog
    LogDebug("HcalSim") << "HcalQie::getEnergy: " << i << " code " << code[i] 
                        << " PE " << work[i];
#endif
  }

  double tmp;
  if (preSamples == 0) {
    tmp = (sum - signalBuckets*baseline) * rescale_ * qToPE * eDepPerPE;
  } else {
    tmp = sum * rescale_ * qToPE;
  }
#ifdef DebugLog
  LogDebug("HcalSim") << "HcalQie::getEnergy: PE " << sum*qToPE << " Energy " 
                      << tmp << " GeV";
#endif
  return tmp;
}

double HcalQie::getShape ( double time ) [private]

Definition at line 288 of file HcalQie.cc.

References gen::k, shape_, and tmp.

Referenced by getCode().

                                    {

  double tmp = 0;
  int    k = (int) (time + 0.5);
  if (k>=0 && k<((int)(shape_.size())-1))
    tmp = 0.5*(shape_[k]+shape_[k+1]);
  
  return tmp;
}

std::vector< double > HcalQie::shape ( void ) [private]

Definition at line 77 of file HcalQie.cc.

References funct::exp(), i, j, gen::k, LogDebug, tmax, and tmp.

Referenced by HcalQie().

                                 {

  // pulse shape time constants in ns
  const float ts1  = 8.;          // scintillation time constants : 1,2,3
  const float ts2  = 10.;           
  const float ts3  = 29.3;         
  const float thpd = 4.;          // HPD current collection drift time
  const float tpre = 5.;          // preamp time constant
  
  const float wd1 = 2.;           // relative weights of decay exponents 
  const float wd2 = 0.7;
  const float wd3 = 1.;

  // HPD starts at I and rises to 2I in thpd of time
  double norm=0.0;
  int j, hpd_siz = (int)(thpd);
  std::vector<double> hpd_drift(hpd_siz);
  for (j=0; j<hpd_siz; j++) {
    double tmp = (double)j + 0.5;
    hpd_drift[j] = 1.0 + tmp/thpd;
    norm += hpd_drift[j];
  }
  // normalize integrated current to 1.0
  for (j=0; j<hpd_siz; j++) {
    hpd_drift[j] /= norm;
  }
  
  // Binkley shape over 6 time constants
  int  preamp_siz=(int)(6*tpre);
  std::vector<double> preamp(preamp_siz);
  norm = 0;
  for (j=0; j<preamp_siz; j++) {
    double tmp = (double)j + 0.5;
    preamp[j]  = tmp*exp(-(tmp*tmp)/(tpre*tpre));
    norm += preamp[j];
  }
  // normalize pulse area to 1.0
  for (j=0; j<preamp_siz; j++) {
    preamp[j] /= norm;
  }

  // ignore stochastic variation of photoelectron emission
  // <...>
  // effective tile plus wave-length shifter decay time over 4 time constants

  int tmax = 6 * (int)ts3;
  std::vector<double> scnt_decay(tmax);
  norm = 0;
  for (j=0; j<tmax; j++) {
    double tmp = (double)j + 0.5;
    scnt_decay[j] = wd1*exp(-tmp/ts1) + wd2*exp(-tmp/ts2) + wd3*exp(-tmp/ts3);
    norm += scnt_decay[j];
  }
  // normalize pulse area to 1.0
  for (j=0; j<tmax; j++) {
    scnt_decay[j] /= norm;
  }

  int nsiz = tmax + hpd_siz + preamp_siz + 1;
  std::vector<double> pulse(nsiz,0.0);  // zeroing output pulse shape
  norm = 0;
  int    i, k;
  for (i=0; i<tmax; i++) {
    int t1 = i; // and ignore jitter from optical path length
    for (j=0; j<hpd_siz; j++) {
      int t2 = t1 + j;
      for (k=0; k<preamp_siz; k++) {
        int   t3 = t2 + k;
        float tmp = scnt_decay[i]*hpd_drift[j]*preamp[k];
        pulse[t3]+= tmp;
        norm     += tmp;
      }
    }
  }
  
  // normalize for 1 GeV pulse height
  edm::LogInfo("HcalSim") << "HcalQie: Convoluted Shape ============== "
                          << "Normalisation " << norm;
  for (i=0; i<nsiz; i++) {
    pulse[i] /= norm;
#ifdef DebugLog
    LogDebug("HcalSim") << "HcalQie: Pulse[" << std::setw(3) << i << "] " 
                        << std::setw(8) << pulse[i];
#endif
  }
  
  return pulse;
}

std::vector< double > HcalQie::weight	(	int	binofmax,
		int	mode,
		int	npre,
		int	numbucket
	)		`[private]`

Definition at line 230 of file HcalQie.cc.

References i, j, LogDebug, and groupFilesInBlocks::temp.

Referenced by HcalQie().

                                                {

  std::vector<double> temp(bucket,0);
  int i;
  for (i=binOfMax-1; i<binOfMax+mode-1; i++)
    temp[i] = 1.;
  if (npre>0) {
    for (i=0; i<npre; i++) {
      int j   = binOfMax-2-i;
      temp[j] =-(double)mode/(double)npre;
    }
  }

#ifdef DebugLog
  int siz = temp.size();
  LogDebug("HcalSim") << "HcalQie: Weights in array of size " << siz 
                      << " and Npre " << npre;
  for (i=0; i<siz; i++)
    LogDebug("HcalSim") << "HcalQie: [Weight[" << i << "] = " << temp[i];
#endif
  return temp;
}