HGCAL ToA for overlapping events

DataFormats/HGCDigi/interface/HGCSample.h

@@ -20,16 +20,17 @@ class HGCSample {
   /**
      @short CTOR
    */
- HGCSample() : value_(0) { }
- HGCSample(uint32_t value) : value_(value) { }
- HGCSample( const HGCSample& o ) : value_(o.value_) { }
+ HGCSample() : value_(0), toaFired_(false) { }
+ HGCSample(uint32_t value) : value_(value), toaFired_(false) { }
+ HGCSample( const HGCSample& o ) : value_(o.value_), toaFired_(o.toaFired_) { }
 
   /**
      @short setters
    */
   void setThreshold(bool thr)           { setWord(thr,  kThreshMask, kThreshShift); }
   void setMode(bool mode)               { setWord(mode, kModeMask,   kModeShift);   }
   void setToA(uint16_t toa)             { setWord(toa,  kToAMask,    kToAShift);    }
+  void setToAValid(bool toaFired)       { toaFired_ = toaFired;  }
   void setData(uint16_t data)           { setWord(data, kDataMask,   kDataShift);   }
   void set(bool thr, bool mode,uint16_t toa, uint16_t data) 
   { 
@@ -56,6 +57,7 @@ class HGCSample {
   uint32_t raw()  const      { return value_;                   }
   bool threshold() const     { return ( (value_ >> kThreshShift) & kThreshMask );  }
   bool mode() const          { return ( (value_ >> kModeShift)   & kModeMask   );  }
+  bool getToAValid() const   { return toaFired_;                }
   uint32_t toa()  const      { return ( (value_ >> kToAShift)    & kToAMask    ); }
   uint32_t data()  const     { return ( (value_ >> kDataShift)   & kDataMask   ); }
   uint32_t operator()()      { return value_;                   }
@@ -77,6 +79,7 @@ class HGCSample {
 
   // a 32-bit word
   uint32_t value_;
+  bool toaFired_;
 };
 
 

RecoLocalCalo/HGCalRecAlgos/interface/HGCalUncalibRecHitRecWeightsAlgo.h

@@ -53,7 +53,7 @@ template<class C> class HGCalUncalibRecHitRecWeightsAlgo
   virtual HGCUncalibratedRecHit makeRecHit( const C& dataFrame ) {
     double amplitude_(-1.),  pedestal_(-1.), jitter_(-1.), chi2_(-1.);
     uint32_t flag = 0;
-    
+
     constexpr int iSample=2; //only in-time sample
     const auto& sample = dataFrame.sample(iSample);
 
@@ -69,16 +69,19 @@ template<class C> class HGCalUncalibRecHitRecWeightsAlgo
         // LG (11/04/2016):
         // offset the TDC upwards to reflect the bin center
 	amplitude_ = ( std::floor(tdcOnsetfC_/adcLSB_) + 1.0 )* adcLSB_ + ( double(sample.data()) + 0.5) * tdcLSB_;
-	jitter_    = double(sample.toa()) * toaLSBToNS_;
+	if(sample.getToAValid()) jitter_    = double(sample.toa()) * toaLSBToNS_;
 	LogDebug("HGCUncalibratedRecHit") << "TDC+: set the charge to: " << amplitude_ << ' ' << sample.data() 
                                           << ' ' << tdcLSB_ << std::endl
                                           << "TDC+: set the ToA to: " << jitter_ << ' ' 
                                           << sample.toa() << ' ' << toaLSBToNS_ << ' '
                                           << " flag=" << flag << std::endl;
       } else {
 	amplitude_ = double(sample.data()) * adcLSB_;
+	if(sample.getToAValid()) jitter_    = double(sample.toa()) * toaLSBToNS_;
 	LogDebug("HGCUncalibratedRecHit") << "ADC+: set the charge to: " << amplitude_ << ' ' << sample.data() 
-                                          << ' ' << adcLSB_ << ' ' << std::endl;
+                                          << ' ' << adcLSB_ << ' ' 
+					  << "TDC+: set the ToA to: " << jitter_ << ' '
+                                          << sample.toa() << ' ' << toaLSBToNS_ << ' '<< std::endl;
       }
     } else {
       amplitude_ = double(sample.data()) * adcLSB_;
@@ -90,7 +93,7 @@ template<class C> class HGCalUncalibRecHitRecWeightsAlgo
     if( ddd_ != nullptr ) {
       HGCalDetId hid(dataFrame.id());
       thickness = ddd_->waferTypeL(hid.wafer());
-    }
+    }    
     amplitude_ = amplitude_/fCPerMIP_[thickness-1];
 
     LogDebug("HGCUncalibratedRecHit") << "Final uncalibrated amplitude : " << amplitude_ << std::endl;

RecoLocalCalo/HGCalRecProducers/plugins/HGCalUncalibRecHitWorkerWeights.cc

@@ -44,7 +44,7 @@ void configureIt(const edm::ParameterSet& conf,
     maker.set_TDCLSB(-1.);
     maker.set_tdcOnsetfC(-1.);
   } 
-    
+
   if( conf.exists(toaLSB_ns) ) {
     maker.set_toaLSBToNS(conf.getParameter<double>(toaLSB_ns));
   } else {

SimCalorimetry/HGCalSimProducers/interface/HGCDigitizer.h

@@ -116,6 +116,8 @@ private :
   uint32_t nEvents_;
 
   std::vector<float> cce_;
+
+  std::map< uint32_t, std::vector< std::pair<float, float> > > hitRefs_bx0;
 };
 
 

SimCalorimetry/HGCalSimProducers/interface/HGCDigitizerBase.h

@@ -46,6 +46,7 @@ class HGCDigitizerBase {
   float keV2fC() const { return keV2fC_; }
   bool toaModeByEnergy() const { return (myFEelectronics_->toaMode()==HGCFEElectronics<DFr>::WEIGHTEDBYE); }
   float tdcOnset() const { return myFEelectronics_->getTDCOnset(); }
+  std::array<float,3> tdcForToaOnset() const { return myFEelectronics_->getTDCForToaOnset(); }
 
   /**
      @short a trivial digitization: sum energies and digitize without noise

SimCalorimetry/HGCalSimProducers/interface/HGCFEElectronics.h

@@ -43,13 +43,27 @@ class HGCFEElectronics
       }
   }
 
+
+  void SetNoiseValues(std::vector<float> noise_fC){
+    for( auto noise : noise_fC ) { noise_fC_.push_back( noise ); }
+  };
+
+  float getTimeJitter(float totalCharge, int thickness){
+    float A2 = jitterNoise2_ns_[thickness-1];
+    float C2 = jitterConstant2_ns_[thickness-1];
+    float X2 = pow((totalCharge/noise_fC_[thickness-1]), 2.);
+    float jitter2 = A2 / X2 + C2;
+    return sqrt(jitter2);
+  };
+
   /**
      @short returns the LSB in MIP currently configured
   */
   float getADClsb()       { return adcLSB_fC_;       }
   float getTDClsb()       { return tdcLSB_fC_;       }  
   float getADCThreshold() { return adcThreshold_fC_; }
   float getTDCOnset()     { return tdcOnset_fC_;     }
+  std::array<float,3> getTDCForToaOnset()    { return tdcForToaOnset_fC_; }
   void setADClsb(float newLSB) { adcLSB_fC_=newLSB; }
 
   /**
@@ -83,13 +97,16 @@ class HGCFEElectronics
   //private members
   uint32_t fwVersion_;
   std::array<float,6> adcPulse_, pulseAvgT_;
+  std::array<float,3> tdcForToaOnset_fC_;
   std::vector<float> tdcChargeDrainParameterisation_;
   float adcSaturation_fC_, adcLSB_fC_, tdcLSB_fC_, tdcSaturation_fC_,
-    adcThreshold_fC_, tdcOnset_fC_, toaLSB_ns_, tdcResolutionInNs_; 
+    adcThreshold_fC_, tdcOnset_fC_, toaLSB_ns_, tdcResolutionInNs_;
+  std::array<float,3> jitterNoise2_ns_, jitterConstant2_ns_;
+  std::vector<float> noise_fC_;
   uint32_t toaMode_;
   bool thresholdFollowsMIP_;
   //caches
-  std::array<bool,hgc::nSamples>  busyFlags, totFlags;
+  std::array<bool,hgc::nSamples>  busyFlags, totFlags, toaFlags;
   hgc::HGCSimHitData newCharge, toaFromToT;
 };
 

SimCalorimetry/HGCalSimProducers/python/hgcalDigitizer_cfi.py

@@ -39,6 +39,10 @@
             adcSaturation_fC  = cms.double(100),
             # the tdc resolution smearing (in picoseconds)
             tdcResolutionInPs = cms.double( 0.001 ),
+            # jitter for timing noise term ns
+            jitterNoise_ns = cms.vdouble(0., 0., 0.),
+            # jitter for timing noise term ns
+            jitterConstant_ns = cms.vdouble(0.00, 0.00, 0.00),
             # LSB for TDC, assuming 12 bit dynamic range to 10 pC
             tdcNbits          = cms.uint32(12),
             # TDC saturation
@@ -48,7 +52,9 @@
             adcThreshold_fC   = cms.double(0.672),
             thresholdFollowsMIP        = cms.bool(thresholdTracksMIP),
             # raise usage of TDC and mode flag (from J. Kaplon)
-            tdcOnset_fC       = cms.double(60) ,
+            tdcOnset_fC       = cms.double(60),
+            # raise usage of TDC for TOA only
+            tdcForToaOnset_fC = cms.vdouble(60., 60., 60.),
             # LSB for time of arrival estimate from TDC in ns
             toaLSB_ns         = cms.double(0.005),
             #toa computation mode (0=by weighted energy, 1=simple threshold)
@@ -92,6 +98,10 @@
             adcSaturation_fC  = cms.double(100),
             # the tdc resolution smearing (in picoseconds)
             tdcResolutionInPs = cms.double( 0.001 ),
+            # jitter for timing noise term ns
+            jitterNoise_ns = cms.vdouble(0., 0., 0.),
+            # jitter for timing noise term ns
+            jitterConstant_ns = cms.vdouble(0.00, 0.00, 0.00),
             # LSB for TDC, assuming 12 bit dynamic range to 10 pC
             tdcNbits          = cms.uint32(12),
             # TDC saturation
@@ -101,7 +111,9 @@
             adcThreshold_fC   = cms.double(0.672),
             thresholdFollowsMIP        = cms.bool(thresholdTracksMIP),
             # raise usage of TDC and mode flag (from J. Kaplon)
-            tdcOnset_fC       = cms.double(60) ,
+            tdcOnset_fC       = cms.double(60), 
+            # raise usage of TDC for TOA only                                                                                                                            
+            tdcForToaOnset_fC = cms.vdouble(60., 60., 60.),
             # LSB for time of arrival estimate from TDC in ns
             toaLSB_ns         = cms.double(0.005),
             #toa computation mode (0=by weighted energy, 1=simple threshold)

SimCalorimetry/HGCalSimProducers/src/HGCDigitizer.cc

@@ -28,6 +28,10 @@ namespace {
 
   constexpr std::array<double,3> occupancyGuesses = { { 0.5,0.2,0.2 } };
 
+  bool comparePairs(const std::pair<float, float>& i, const std::pair<float, float>& j){
+    return i.second < j.second;
+  }
+
   float getPositionDistance(const HGCalGeometry* geom, const DetId& id) {
     return geom->getPosition(id).mag();
   }
@@ -36,6 +40,20 @@ namespace {
     return geom->getGeometry(id)->getPosition().mag();
   }
 
+  int getCellThickness(const HGCalGeometry* geom, const DetId& detid ) {
+    const auto& topo     = geom->topology();
+    const auto& dddConst = topo.dddConstants();
+    uint32_t id(detid.rawId());
+    HGCalDetId hid(id);
+    int wafer = HGCalDetId(id).wafer();
+    int waferTypeL = dddConst.waferTypeL(wafer);
+    return waferTypeL;
+  }
+
+  int getCellThickness(const HcalGeometry* geom, const DetId& detid ) {
+    return 1;
+  }
+
   void getValidDetIds(const HGCalGeometry* geom, std::unordered_set<DetId>& valid) {
     const std::vector<DetId>& ids = geom->getValidDetIds();
     valid.reserve(ids.size());
@@ -187,6 +205,7 @@ void HGCDigitizer::initializeEvent(edm::Event const& e, edm::EventSetup const& e
 //
 void HGCDigitizer::finalizeEvent(edm::Event& e, edm::EventSetup const& es, CLHEP::HepRandomEngine* hre)
 {
+  hitRefs_bx0.clear();
 
   const CaloSubdetectorGeometry* theGeom = ( nullptr == gHGCal_ ? 
 					     static_cast<const CaloSubdetectorGeometry*>(gHcal_) : 
@@ -296,21 +315,22 @@ void HGCDigitizer::accumulate(edm::Handle<edm::PCaloHitContainer> const &hits,
 
   //configuration to apply for the computation of time-of-flight
   bool weightToAbyEnergy(false);
-  float tdcOnset(0.f),keV2fC(0.f);
+  std::array<float, 3> tdcForToaOnset{ {0.f, 0.f, 0.f} };
+  float keV2fC(0.f);
   switch( mySubDet_ ) {
   case ForwardSubdetector::HGCEE:
     weightToAbyEnergy = theHGCEEDigitizer_->toaModeByEnergy();
-    tdcOnset          = theHGCEEDigitizer_->tdcOnset();
+    tdcForToaOnset    = theHGCEEDigitizer_->tdcForToaOnset();
     keV2fC            = theHGCEEDigitizer_->keV2fC();
     break;
   case ForwardSubdetector::HGCHEF:
     weightToAbyEnergy = theHGCHEfrontDigitizer_->toaModeByEnergy();
-    tdcOnset          = theHGCHEfrontDigitizer_->tdcOnset();
+    tdcForToaOnset    = theHGCHEfrontDigitizer_->tdcForToaOnset();
     keV2fC            = theHGCHEfrontDigitizer_->keV2fC();
     break;
   case ForwardSubdetector::HGCHEB:
     weightToAbyEnergy = theHGCHEbackDigitizer_->toaModeByEnergy();
-    tdcOnset          = theHGCHEbackDigitizer_->tdcOnset();
+    tdcForToaOnset    = theHGCHEbackDigitizer_->tdcForToaOnset();
     keV2fC            = theHGCHEbackDigitizer_->keV2fC();     
     break;
   default:
@@ -363,7 +383,7 @@ void HGCDigitizer::accumulate(edm::Handle<edm::PCaloHitContainer> const &hits,
     //accumulate in 15 buckets of 25ns (9 pre-samples, 1 in-time, 5 post-samples)
     const float tof = toa-dist2center/refSpeed_+tofDelay_ ;
     const int itime= std::floor( tof/bxTime_ ) + 9;
-      
+
     //no need to add bx crossing - tof comes already corrected from the mixing module
     //itime += bxCrossing;
     //itime += 9;
@@ -375,32 +395,59 @@ void HGCDigitizer::accumulate(edm::Handle<edm::PCaloHitContainer> const &hits,
 
     (simHitIt->second).hit_info[0][itime] += charge;
     float accCharge=(simHitIt->second).hit_info[0][itime];
-
+
+
+    //working version with pileup only for in-time hits
+    int waferThickness = getCellThickness(geom,id);
+    float accChargeForToA = 0.f;
+    bool orderChanged = false;
+    if(itime == 9){
+      if(hitRefs_bx0[id].empty()){
+	hitRefs_bx0[id].push_back(std::pair<float, float>(charge, tof));
+	accChargeForToA += charge;
+      }
+      else if(tof <= hitRefs_bx0[id].back().second){
+	hitRefs_bx0[id].push_back(std::pair<float, float>(charge, tof));
+	std::sort(hitRefs_bx0[id].begin(), hitRefs_bx0[id].end(), comparePairs);
+	for(const auto& step : hitRefs_bx0[id]){
+	  accChargeForToA += step.first;
+	  if(accChargeForToA > tdcForToaOnset[waferThickness-1] && step.second != hitRefs_bx0[id].back().second){
+	    while(step != hitRefs_bx0[id].back())  hitRefs_bx0[id].pop_back();
+	    break;
+	  }
+	}
+	orderChanged = true;
+      }
+      else{
+        if(accCharge - charge <= tdcForToaOnset[waferThickness-1]){
+          hitRefs_bx0[id].push_back(std::pair<float, float>(charge, tof));
+          accChargeForToA = accCharge;
+        }
+      }
+    }
+
     //time-of-arrival (check how to be used)
     if(weightToAbyEnergy) (simHitIt->second).hit_info[1][itime] += charge*tof;
-    else if((simHitIt->second).hit_info[1][itime]==0) {	
-      if( accCharge>tdcOnset)
-	{
-	  //extrapolate linear using previous simhit if it concerns to the same DetId
-	  float fireTDC=tof;
-	  if(i>0)
-	    {
-	      uint32_t prev_id  = std::get<1>(hitRefs[i-1]);
-	      if(prev_id==id)
-		{
-		  float prev_toa    = std::get<2>(hitRefs[i-1]);
-		  float prev_tof(prev_toa-dist2center/refSpeed_+tofDelay_);
-		  //float prev_charge = std::get<3>(hitRefs[i-1]);
-		  float deltaQ2TDCOnset = tdcOnset-((simHitIt->second).hit_info[0][itime]-charge);
-		  float deltaQ          = charge;
-		  float deltaT          = (tof-prev_tof);
-		  fireTDC               = deltaT*(deltaQ2TDCOnset/deltaQ)+prev_tof;
-		}		  
-	    }
-
-	  (simHitIt->second).hit_info[1][itime]=fireTDC;
+    else if(accChargeForToA > tdcForToaOnset[waferThickness-1] &&
+	    ((simHitIt->second).hit_info[1][itime] == 0 || orderChanged == true) ){
+      float fireTDC = hitRefs_bx0[id].back().second;
+      if (hitRefs_bx0[id].size() > 1){
+	float chargeBeforeThr = 0.f;
+	float tofchargeBeforeThr = 0.f;
+	for(const auto& step : hitRefs_bx0[id]){
+	  if(step.first + chargeBeforeThr <= tdcForToaOnset[waferThickness-1]){
+	    chargeBeforeThr += step.first;
+	    tofchargeBeforeThr = step.second;
+	  }
+	  else break;
 	}
+	float deltaQ = accChargeForToA - chargeBeforeThr;
+	float deltaTOF = fireTDC - tofchargeBeforeThr;
+	fireTDC = (tdcForToaOnset[waferThickness-1] - chargeBeforeThr) * deltaTOF / deltaQ + tofchargeBeforeThr;
+      }
+      (simHitIt->second).hit_info[1][itime] = fireTDC;                                                                  
     }
+
   }
   hitRefs.clear();
 }

SimCalorimetry/HGCalSimProducers/src/HGCDigitizerBase.cc

@@ -73,6 +73,7 @@ HGCDigitizerBase<DFr>::HGCDigitizerBase(const edm::ParameterSet& ps) {
   }
   edm::ParameterSet feCfg = myCfg_.getParameter<edm::ParameterSet>("feCfg");
   myFEelectronics_        = std::unique_ptr<HGCFEElectronics<DFr> >( new HGCFEElectronics<DFr>(feCfg) );
+  myFEelectronics_->SetNoiseValues(noise_fC_); 
 }
 
 template<class DFr>
@@ -120,6 +121,7 @@ void HGCDigitizerBase<DFr>::runSimple(std::unique_ptr<HGCDigitizerBase::DColl> &
 
       //add noise (in fC)
       //we assume it's randomly distributed and won't impact ToA measurement
+      //also assume that it is related to the charge path only and that noise fluctuation for ToA circuit be handled separately
       totalCharge += std::max( (float)CLHEP::RandGaussQ::shoot(engine,0.0,cell.size*noise_fC_[cell.thickness-1]) , 0.f );
       if(totalCharge<0.f) totalCharge=0.f;