TDC changes on top of CMSSW_7_0_X

CalibCalorimetry/HcalAlgos/src/HcalPulseShapes.cc

@@ -286,19 +286,145 @@ void HcalPulseShapes::computeHFShape() {
 
 void HcalPulseShapes::computeSiPMShape()
 {
-  unsigned int nbin = 100; // to avoid big drop of integral for previous 512
-                       // due to negative afterpulse (May 6, 2013. S.Abdullin) 
+
+  unsigned int nbin = 128; 
+
+//From Jake Anderson: numberical convolution of SiPMs  WLC shapes
+  std::vector<float> nt = {
+    2.782980485851731e-6,
+    4.518134885954626e-5,
+    2.7689305197392056e-4,
+    9.18328418900969e-4,
+    .002110072599166349,
+    .003867856860331454,
+    .006120046224897771,
+    .008754774090536956,
+    0.0116469503358586,
+    .01467007449455966,
+    .01770489955229477,
+    .02064621450689512,
+    .02340678093764222,
+    .02591874610854916,
+    .02813325527435303,
+    0.0300189241965647,
+    .03155968107671164,
+    .03275234052577155,
+    .03360415306318798,
+    .03413048377960748,
+    .03435270899678218,
+    .03429637464659661,
+    .03398962975487166,
+    .03346192884394954,
+    .03274298516247742,
+    .03186195009136525,
+    .03084679116113031,
+    0.0297238406141036,
+    .02851748748929785,
+    .02724998816332392,
+    .02594137274487424,
+    .02460942736731527,
+    .02326973510736116,
+    .02193576080366117,
+    0.0206189674254987,
+    .01932895378564653,
+    0.0180736052958666,
+    .01685925112650875,
+    0.0156908225633535,
+    .01457200857138456,
+    .01350540559602467,
+    .01249265947824805,
+    .01153459805300423,
+    .01063135355597282,
+    .009782474412011936,
+    .008987026319784546,
+    0.00824368281357106,
+    .007550805679909604,
+    .006906515742762193,
+    .006308754629755056,
+    .005755338185695127,
+    .005244002229973356,
+    .004772441359900532,
+    .004338341490928299,
+    .003939406800854143,
+    0.00357338171220501,
+    0.0032380685079891,
+    .002931341133259233,
+    .002651155690306086,
+    .002395558090237333,
+    .002162689279320922,
+    .001950788415487319,
+    .001758194329648101,
+    .001583345567913682,
+    .001424779275191974,
+    .001281129147671334,
+    0.00115112265163774,
+    .001033577678808199,
+    9.273987838127585e-4,
+    8.315731274976846e-4,
+    7.451662302008696e-4,
+    6.673176219006913e-4,
+    5.972364609644049e-4,
+    5.341971801529036e-4,
+    4.775352065178378e-4,
+    4.266427928961177e-4,
+    3.8096498904225923e-4,
+    3.3999577417327287e-4,
+    3.032743659102713e-4,
+    2.703817158798329e-4,
+    2.4093719775272793e-4,
+    2.145954900503894e-4,
+    1.9104365317752797e-4,
+    1.6999839784346724e-4,
+    1.5120354022478893e-4,
+    1.3442763782650755e-4,
+    1.1946179895521507e-4,
+    1.0611765796993575e-4,
+    9.422550797617687e-5,
+    8.363258233342666e-5,
+    7.420147621931836e-5,
+    6.580869950304933e-5,
+    5.834335229919868e-5,
+    5.17059147771959e-5,
+    4.5807143072062634e-5,
+    4.0567063461299446e-5,
+    3.591405732740723e-5,
+    3.178402980354131e-5,
+    2.811965539165646e-5,
+    2.4869694240316126e-5,
+    2.1988373166730962e-5,
+    1.9434825899529382e-5,
+    1.717258740121378e-5,
+    1.5169137499243157e-5,
+    1.339548941011129e-5,
+    1.1825819079078403e-5,
+    1.0437131581057595e-5,
+    9.208961130078894e-6,
+    8.12310153137994e-6,
+    7.163364176588591e-6,
+    6.315360932244386e-6,
+    5.566309502463164e-6,
+    4.904859063429651e-6,
+    4.320934164082596e-6,
+    3.8055950719111903e-6,
+    3.350912911083174e-6,
+    2.9498580949517117e-6,
+    2.596200697612328e-6,
+    2.2844215378879293e-6,
+    2.0096328693141094e-6,
+    1.7675076766686654e-6,
+    1.5542166787225756e-6,
+    1.366372225473431e-6,
+    1.200978365778838e-6,
+    1.0553864128982371e-6,
+    9.272554464808518e-7,
+    8.145171945902259e-7,
+    7.153448381918271e-7
+  };
+
   siPMShape_.setNBin(nbin);
-  std::vector<float> nt(nbin,0.0);  //
 
   double norm = 0.;
   for (unsigned int j = 0; j < nbin; ++j) {
-    if (j <= 31.)
-      nt[j] = 2.15*j;
-    else if ((j > 31) && (j <= 96))
-      nt[j] = 102.1 - 1.12*j;
-    else
-      nt[j] = 0.0076*j - 6.4;
     norm += (nt[j]>0) ? nt[j] : 0.;
   }
 
@@ -308,6 +434,11 @@ void HcalPulseShapes::computeSiPMShape()
   }
 }
 
+double HcalPulseShapes::gexp(double t, double A, double c, double t0, double s) {
+  static double const root2(sqrt(2));
+  return -A*0.5*exp(c*t+0.5*c*c*s*s-c*s)*(erf(-0.5*root2/s*(t-t0+c*s*s))-1);
+}
+
 
 const HcalPulseShapes::Shape &
 HcalPulseShapes::getShape(int shapeType) const

CalibFormats/CaloObjects/src/CaloSamples.cc

@@ -102,20 +102,20 @@ CaloSamples::setBlank() // keep id, presamples, size but zero out data
 }
 
 std::ostream& operator<<(std::ostream& s, const CaloSamples& samples) {
-  s << "DetId=" << samples.id();
+  s << "DetId " << samples.id();
   // print out every so many precise samples
   float preciseStep = samples.preciseSize()/samples.size();
   s << ", "<<  samples.size() << " samples";
   if (preciseStep > 0) 
     s << ", " << samples.preciseSize() << " preciseSamples" 
-      << ", " << preciseStep << " precise step";
+      << ", " << preciseStep << " preciseStep";
   s << '\n';
   for (int i=0; i<samples.size(); i++) {
-    s << i << ":" << samples[i];
-    int precisei = i*preciseStep;
-    if(precisei < samples.preciseSize()) {
-      s << " " << samples.preciseAt(precisei) ;
-    }
+    s << i << ":" << samples[i] << " precise:";
+    int precise_start(i*preciseStep), precise_end(precise_start + preciseStep);
+    for (int j(precise_start); ((j<precise_end) && (j< samples.preciseSize()));
+	 ++j)
+      s << " " << samples.preciseAt(j);
     s << std::endl;
   }
   return s;

RecoLocalCalo/HcalRecAlgos/src/HcalTDCReco.cc

@@ -1,48 +1,57 @@
-#include "RecoLocalCalo/HcalRecAlgos/src/HcalTDCReco.h"
-#include "DataFormats/HcalDigi/interface/HcalUpgradeDataFrame.h"
-#include "DataFormats/HcalRecHit/interface/HBHERecHit.h"
-
-HcalTDCReco::HcalTDCReco() 
-{}
-
-void HcalTDCReco::reconstruct(const HcalUpgradeDataFrame & digi, 
-                              HBHERecHit & recHit) const
-{
-  int n = digi.size();
-  double risingTime = -999.;
-  double fallingTime = -999.;
-  int signalBX = 5;              // NB: HARDWIRED !!!
-  int nbins = 50; // as in HcalTDCParameters.h (SimCalorimetry/HcalSimAlgos)
-  // start at bunch crossing 3 by default
-  for(int i = 3; i < n; ++i)
-  {
-    unsigned tdc = digi.tdc(i);
-
-    /* 
-    unsigned rising = tdc & 0x7F;
-    unsigned falling = (tdc >> 7) & 0x7F;
-    */
-    // temporary "unpacking" instead, which directly corresponds to 
-    // SimCalorimetry/HcalSimAlgos/src/HcalTDC.cc  (nibs = 50...)
-    // packedTDC = TDC_RisingEdge + (tdcBins*2) * TDC_FallingEdge;
-    unsigned rising  =  tdc%100; 
-    unsigned falling =  tdc/100;
-
-    // only set the first time, avoiding "special" codes 
-    if(risingTime < -998. && rising != 62 && rising != 63) {
-      risingTime = rising*25./nbins + (i-signalBX)*25.;
-    }
-    if(fallingTime < -998. && falling != 62 && falling != 63) {
-      fallingTime = falling*25./nbins + (i-signalBX)*25.;
-    }    
-
-    /*
-    std::cout << " digi.tdc[" << i << "] = " << tdc 
-	      << "  rising = " << rising << "  falling = " << falling
-	      << "   Rt = " << risingTime
-	      << "   Ft = " << fallingTime
-	      << std::endl;
-    */
-  }
-  recHit = HBHERecHit(recHit.id(), recHit.energy(), risingTime, fallingTime);
-}
+#include "RecoLocalCalo/HcalRecAlgos/src/HcalTDCReco.h"
+#include "DataFormats/HcalDigi/interface/HcalUpgradeDataFrame.h"
+#include "DataFormats/HcalRecHit/interface/HBHERecHit.h"
+
+HcalTDCReco::HcalTDCReco() 
+{}
+
+void HcalTDCReco::reconstruct(const HcalUpgradeDataFrame & digi, 
+                              HBHERecHit & recHit) const
+{
+  int n = digi.size();
+  double risingTime = -999.;
+  double fallingTime = -999.;
+  int signalBX = 4;              // NB: HARDWIRED !!!
+  int nbins = 50; // as in HcalTDCParameters.h (SimCalorimetry/HcalSimAlgos)
+  int direction(-1), stepSize(1); // where to go after signalBX
+  // start at bunch crossing 3 by default
+  int i(signalBX);
+
+  // for(int i = 3; i < n; ++i)
+  while ((i > 2) && (i < 8) && (i < n) && 
+	 ((risingTime < -998.) || (fallingTime < 998.)))
+  {
+    unsigned tdc = digi.tdc(i);
+
+    /* 
+    unsigned rising = tdc & 0x7F;
+    unsigned falling = (tdc >> 7) & 0x7F;
+    */
+    // temporary "unpacking" instead, which directly corresponds to 
+    // SimCalorimetry/HcalSimAlgos/src/HcalTDC.cc  (nibs = 50...)
+    // packedTDC = TDC_RisingEdge + (tdcBins*2) * TDC_FallingEdge;
+    unsigned rising  =  tdc%100; 
+    unsigned falling =  tdc/100;
+
+    // only set the first time, avoiding "special" codes 
+    if(risingTime < -998. && rising != 62 && rising != 63) {
+      risingTime = rising*25./nbins + (i-signalBX)*25.;
+    }
+    if(((fallingTime < -998.) || (fallingTime < risingTime)) && 
+       (falling != 62) && (falling != 63)) {
+      fallingTime = falling*25./nbins + (i-signalBX)*25.;
+    }  
+
+    i += direction*stepSize;
+    ++stepSize;
+    direction *= -1;
+    /*
+    std::cout << " digi.tdc[" << i << "] = " << tdc 
+	      << "  rising = " << rising << "  falling = " << falling
+	      << "   Rt = " << risingTime
+	      << "   Ft = " << fallingTime
+	      << std::endl;
+    */
+  }
+  recHit = HBHERecHit(recHit.id(), recHit.energy(), risingTime, fallingTime);
+}

RecoLocalCalo/HcalRecProducers/python/HBHEUpgradeReconstructor_cfi.py

@@ -1,14 +1,14 @@
-import FWCore.ParameterSet.Config as cms
-
-hbheUpgradeReco = cms.EDProducer("HcalSimpleReconstructor",
-    correctionPhaseNS = cms.double(-5.0),  
-    digiLabel = cms.InputTag("simHcalUnsuppressedDigis","HBHEUpgradeDigiCollection"),
-    Subdetector = cms.string('upgradeHBHE'),
-    correctForPhaseContainment = cms.bool(True),
-    correctForTimeslew = cms.bool(True),
-    dropZSmarkedPassed = cms.bool(True),
-    firstSample = cms.int32(4),
-    samplesToAdd = cms.int32(2),
-    tsFromDB = cms.bool(True)
-)
-
+import FWCore.ParameterSet.Config as cms
+
+hbheUpgradeReco = cms.EDProducer("HcalSimpleReconstructor",
+    correctionPhaseNS = cms.double(5.0),  
+    digiLabel = cms.InputTag("simHcalUnsuppressedDigis","HBHEUpgradeDigiCollection"),
+    Subdetector = cms.string('upgradeHBHE'),
+    correctForPhaseContainment = cms.bool(True),
+    correctForTimeslew = cms.bool(False),
+    dropZSmarkedPassed = cms.bool(True),
+    firstSample = cms.int32(4),
+    samplesToAdd = cms.int32(2),
+    tsFromDB = cms.bool(True)
+)
+

SimCalorimetry/HcalSimAlgos/interface/HcalSiPMShape.h

@@ -15,8 +15,15 @@ class HcalSiPMShape : public CaloVShape {
 
   virtual double operator() (double time) const;
 
-    virtual double       timeToRise()         const {return 33.;}
+  virtual double timeToRise() const {return 3.5;}
+
+  static double gexp(double t, double A, double c, double t0, double s);
+  static double gexpIndefIntegral(double t, double A, double c, double t0, 
+				  double s);
+  static double gexpIntegral0Inf(double A, double c, double t0, double s);
+
 protected:
+  virtual double analyticPulseShape(double t) const;
   void computeShape();
 
 private:

SimCalorimetry/HcalSimAlgos/interface/HcalSimParameters.h

@@ -35,6 +35,7 @@ class HcalSimParameters : public CaloSimParameters
   double timeSmearRMS(double ampl) const;
 
   int pixels() const {return thePixels;}
+  bool doSiPMSmearing() const { return theSiPMSmearing; }
 
   friend class HcalSimParameterMap;
 
@@ -45,6 +46,7 @@ class HcalSimParameters : public CaloSimParameters
   std::vector<double> theSamplingFactors;
   std::vector<double> thePE2fCByRing;
   int thePixels;
+  bool theSiPMSmearing;
   bool doTimeSmear_;
   HcalTimeSmearSettings theSmearSettings;
 };

SimCalorimetry/HcalSimAlgos/interface/HcalTDC.h

@@ -1,3 +1,4 @@
+// -*- mode: c++ -*- 
 #ifndef HcalSimAlgos_HcalTDC_h
 #define HcalSimAlgos_HcalTDC_h
 
@@ -26,11 +27,15 @@ class HcalTDC {
 
 private:
   double getThreshold(const HcalGenericDetId & detId) const;
+  double getHysteresisThreshold(double nominal) const;
+
   HcalTDCParameters theTDCParameters;
   const HcalDbService * theDbService;
 
   unsigned int theDAC;
 
+  double const lsb;
+
   mutable CLHEP::RandGaussQ * theRandGaussQ;
 };
 

SimCalorimetry/HcalSimAlgos/src/HcalElectronicsSim.cc

@@ -67,6 +67,7 @@ void HcalElectronicsSim::analogToDigital(CaloSamples & lf, ZDCDataFrame & result
 void HcalElectronicsSim::analogToDigital(CaloSamples & lf, 
 					 HcalUpgradeDataFrame & result) {
   convert<HcalUpgradeDataFrame>(lf, result);
+//   std::cout << HcalDetId(lf.id()) << ' ' << lf;
   theTDC.timing(lf, result);
 }