HB radiation damage model

CalibCalorimetry/HcalPlugins/python/Hcal_Conditions_forGlobalTag_cff.py

@@ -12,8 +12,10 @@
 toGet = cms.untracked.vstring('GainWidths'),
 #--- the following 5 parameters can be omitted in case of regular Geometry 
     iLumi = cms.double(-1.),                 # for Upgrade: fb-1
-    HERecalibration = cms.uint32(0),     # 1 for Upgrade (default aging scenario)
-    HEreCalibCutoff = cms.double(20.),       # if above is True  
+    HERecalibration = cms.uint32(0),         # 1 for Upgrade (default aging scenario)
+    HEreCalibCutoff = cms.double(20.),       # if above is True
+    HBRecalibration = cms.uint32(0),         # 1 for Upgrade (default aging scenario)
+    HBreCalibCutoff = cms.double(20.),       # if above is True
     HFRecalibration = cms.bool(False),       # True for Upgrade
     GainWidthsForTrigPrims = cms.bool(False) # True Upgrade    
 )

CalibCalorimetry/HcalPlugins/src/HBRecalibration.cc

@@ -0,0 +1,115 @@
+///////////////////////////////////////////////////////////////////////////////
+// File: HBRecalibration.cc
+// Description: simple helper class containing tabulater/parameterized 
+//              function for HE damage recovery for Upgrade studies  
+//              evaluated using DataFormats/HcalCalibObjects/HBDarkening   
+///////////////////////////////////////////////////////////////////////////////
+
+#include "HBRecalibration.h"
+
+HBRecalibration::HBRecalibration(double integrated_lumi, double cutoff, unsigned int scenario):
+iLumi(integrated_lumi),cutoff_(cutoff),darkening(scenario)
+{ }
+
+HBRecalibration::~HBRecalibration() { }
+
+void HBRecalibration::setDsegm( const std::vector<std::vector<int>>& m_segmentation) 
+{
+
+  //  std::cout << std::endl << " HBRecalibration->setDsegm" << std::endl;
+
+  for (unsigned int ieta = 0; ieta < HBDarkening::nEtaBins; ieta++) {
+    //    std::cout << "["<< ieta << "]  ieta =" << ieta + 16 << "  ";
+    for(unsigned int ilay = 0; ilay < HBDarkening::nScintLayers; ilay++) {
+      dsegm[ieta][ilay] = m_segmentation[ieta+1][ilay]; // 0 not used
+      //      std::cout << dsegm [ieta][ilay];
+    }
+    //    std::cout << std::endl;
+  }
+
+  initialize();
+
+}
+
+double HBRecalibration::getCorr(int ieta, int idepth) {
+
+  //  int init_ieta = ieta;
+  ieta = abs(ieta)-1; // 0 - 15
+
+  if(corr[ieta][idepth] > cutoff_) return cutoff_;
+  else return corr[ieta][idepth];
+}
+
+
+void HBRecalibration::initialize() {
+
+  double dval[HBDarkening::nEtaBins][nDepths];  // conversion of lval into depths-averaged values - denominator (including degradation for iLumi) 
+  double nval[HBDarkening::nEtaBins][nDepths];  // conversion of lval into depths-averaged values - numerator (no degradation)
+
+  for (unsigned int j = 0; j < HBDarkening::nEtaBins; j++) {
+    for (unsigned int k = 0; k < nDepths; k++) {
+      dval[j][k] = 0.0;
+      nval[j][k] = 0.0;
+    }
+  }
+
+  double lval[HBDarkening::nEtaBins][HBDarkening::nScintLayers]    // raw table of mean energy in each layer for each ieta at 0 lumi
+    = {
+      {2.362808,1.575159,1.283007,1.026073,0.834189,0.702393,0.566008,0.484473,0.402106,0.306254,0.251159,0.199382,0.156932,0.132067,0.099506,0.080853,0.118480}, //tower 1
+      {2.397443,1.616431,1.301160,1.078375,0.882232,0.705615,0.577277,0.472366,0.383500,0.326740,0.265406,0.208601,0.169150,0.124831,0.103368,0.078420,0.117037}, //tower 2
+      {2.475831,1.654173,1.322021,1.033156,0.850223,0.684293,0.551837,0.467637,0.377493,0.329712,0.254339,0.203073,0.165079,0.123208,0.102514,0.079879,0.112186}, //tower 3
+      {2.462020,1.606079,1.282836,1.010049,0.855179,0.669815,0.556634,0.452061,0.377653,0.311965,0.255425,0.194654,0.159325,0.117479,0.092718,0.069317,0.105891}, //tower 4
+      {2.707523,1.747932,1.373208,1.106487,0.866754,0.728357,0.563662,0.469271,0.374237,0.303459,0.241259,0.187083,0.140809,0.116930,0.081804,0.064648,0.101095}, //tower 5
+      {2.678347,1.741005,1.317389,1.013474,0.823785,0.710703,0.551169,0.429373,0.347780,0.274428,0.222277,0.166950,0.131702,0.101153,0.075034,0.060390,0.092811}, //tower 6
+      {2.809996,1.729876,1.328158,1.050613,0.820268,0.669393,0.516984,0.411677,0.331940,0.272142,0.197038,0.155147,0.127178,0.094173,0.069392,0.054302,0.071747}, //tower 7
+      {2.858155,1.770711,1.355659,1.047950,0.851036,0.657216,0.526521,0.416379,0.319870,0.247920,0.188156,0.145095,0.110619,0.080743,0.063796,0.052655,0.071961}, //tower 8
+      {3.041316,1.877900,1.418290,1.082000,0.840637,0.676057,0.505557,0.402506,0.301596,0.231630,0.182558,0.138136,0.106255,0.073500,0.056109,0.041644,0.045892}, //tower 9
+      {3.142461,1.817359,1.363827,1.013841,0.768494,0.603310,0.463155,0.368469,0.282965,0.218877,0.152383,0.118167,0.079790,0.053056,0.038893,0.031578,0.040494}, //tower 10
+      {3.294945,1.854570,1.367346,1.008908,0.769117,0.594254,0.445583,0.335300,0.248729,0.190494,0.137710,0.099664,0.071549,0.053245,0.037729,0.025820,0.033520}, //tower 11
+      {3.579348,1.951865,1.393643,1.009883,0.745045,0.555898,0.424502,0.313736,0.220994,0.156913,0.109682,0.074274,0.052242,0.039518,0.029820,0.016796,0.028697}, //tower 12
+      {3.752190,2.001947,1.431959,1.057370,0.750016,0.550919,0.410243,0.300922,0.201059,0.151015,0.111034,0.074270,0.049154,0.034618,0.025749,0.018542,0.025501}, //tower 13
+      {4.057676,2.074200,1.423406,1.000701,0.745990,0.536405,0.375057,0.278802,0.192733,0.128711,0.094431,0.062570,0.051008,0.034336,0.026740,0.015083,0.002232}, //tower 14
+      {4.252095,2.160939,1.492043,1.021666,0.740503,0.534830,0.378022,0.276588,0.204727,0.156732,0.111826,0.066944,0.045150,0.031974,0.019658,0.001816,0.000000}, //tower 15
+      {0.311054,0.185613,2.806550,1.810331,0.996787,0.539745,0.285648,0.136337,0.059499,0.007867,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000} //tower 16
+    };
+
+  // coverting energy values from layers into depths 
+
+  //  std::cout << std::endl << " >>> DVAL evaluation " << std::endl;
+
+  for (unsigned int ieta = 0; ieta < HBDarkening::nEtaBins; ieta++) {
+
+    //fill sum(means(layer,0)) and sum(means(layer,lumi)) for each depth
+    for(unsigned int ilay = 0; ilay < HBDarkening::nScintLayers; ilay++) {
+      int idepth = dsegm[ieta][ilay]; // idepth = 0 - not used!
+      nval[ieta][idepth] += lval[ieta][ilay];
+      dval[ieta][idepth] += lval[ieta][ilay]*darkening.degradation(iLumi,ieta+1,ilay); //be careful of eta and layer numbering
+
+      /*
+      std::cout << "ilay " << ilay << " -> idepth " << idepth   
+              << "  + lval[" << ieta << "][" << ilay << "]"
+       << " " << lval[ieta][ilay] <<  std::endl;       
+      */
+    }
+
+    //compute factors, w/ safety checks
+    for(unsigned int idepth = 0; idepth < nDepths; idepth++){
+      if(dval[ieta][idepth] > 0) corr[ieta][idepth] = nval[ieta][idepth]/dval[ieta][idepth];
+      else corr[ieta][idepth] = 1.0;
+
+      if(corr[ieta][idepth] < 1.0) corr[ieta][idepth] = 1.0;
+      /*
+	if (idepth > 0 && idepth <= 3) {
+	std::cout << "nval[" << ieta << "][" << idepth << "]"
+	<< " = " << nval[ieta][idepth] << " - "
+	  	    << "dval["<< ieta << "][" << idepth << "]"
+	  	    << " = " << dval[ieta][idepth] 
+		    << "    corr = " << corr[ieta][idepth] << std::endl;
+		    }
+      */
+    }
+
+  }
+
+
+}

CalibCalorimetry/HcalPlugins/src/HBRecalibration.h

@@ -0,0 +1,41 @@
+#ifndef CalibCalorimetry_HBRecalibration_h
+#define CalibCalorimetry_HBRecalibration_h
+//
+// Simple class with tabulated/parameterized function 
+// to compansate for darkening attenuation of HE scintillators   
+// in Upgrade conditions
+// Evaluated on the basis of DataFormats/HcalCalibObjects/HBDarkening by K.Pedro (Maryland)
+// correction = f (integrated lumi, depth, ieta)
+//
+
+#include <cmath>
+#include <vector>
+#include <iostream>
+#include "DataFormats/HcalCalibObjects/interface/HBDarkening.h"
+
+class HBRecalibration {
+
+public:
+  HBRecalibration(double integrated_lumi, double cutoff, unsigned int scenario);
+  ~HBRecalibration();
+
+  double getCorr(int ieta, int idepth);
+  void  setDsegm(const std::vector<std::vector<int> >& m_segmentation);
+
+private:
+  // max number of HB recalibration depths
+  static const unsigned int nDepths = 7; 
+
+  void initialize();
+  double iLumi;
+  double cutoff_;
+  HBDarkening darkening;
+
+ // Tabulated mean energy values per layer and per depth
+  double dsegm[HBDarkening::nEtaBins][HBDarkening::nScintLayers];
+  double  corr[HBDarkening::nEtaBins][nDepths];
+
+};
+
+
+#endif // HBRecalibration_h

CalibCalorimetry/HcalPlugins/src/HERecalibration.cc

@@ -21,7 +21,7 @@ void HERecalibration::setDsegm( const std::vector<std::vector<int>>& m_segmentat
   for (unsigned int ieta = 0; ieta < HEDarkening::nEtaBins; ieta++) {
     //    std::cout << "["<< ieta << "]  ieta =" << ieta + 16 << "  ";
     for(unsigned int ilay = 0; ilay < HEDarkening::nScintLayers; ilay++) {
-      dsegm[ieta][ilay] = m_segmentation[ieta+15][ilay]; // 0 not used
+      dsegm[ieta][ilay] = m_segmentation[ieta+16][ilay]; // 0 not used
       //      std::cout << dsegm [ieta][ilay];
     }
     //    std::cout << std::endl;

CalibCalorimetry/HcalPlugins/src/HcalHardcodeCalibrations.cc

@@ -108,7 +108,7 @@ namespace {
 
 }
 
-HcalHardcodeCalibrations::HcalHardcodeCalibrations ( const edm::ParameterSet& iConfig ): he_recalibration(0), hf_recalibration(0), setHEdsegm(false)
+HcalHardcodeCalibrations::HcalHardcodeCalibrations ( const edm::ParameterSet& iConfig ): he_recalibration(0), hb_recalibration(0), hf_recalibration(0), setHEdsegm(false), setHBdsegm(false)
 {
   edm::LogInfo("HCAL") << "HcalHardcodeCalibrations::HcalHardcodeCalibrations->...";
 
@@ -124,11 +124,16 @@ HcalHardcodeCalibrations::HcalHardcodeCalibrations ( const edm::ParameterSet& iC
 
   if( iLumi > 0.0 ) {
     unsigned he_recalib = iConfig.getParameter<unsigned>("HERecalibration");
+    unsigned hb_recalib = iConfig.getParameter<unsigned>("HBRecalibration");
     bool hf_recalib = iConfig.getParameter<bool>("HFRecalibration");
     if(he_recalib) {
       double cutoff = iConfig.getParameter<double>("HEreCalibCutoff"); 
       he_recalibration = new HERecalibration(iLumi,cutoff,he_recalib);
     }
+    if(hb_recalib) {
+      double cutoff = iConfig.getParameter<double>("HBreCalibCutoff"); 
+      hb_recalibration = new HBRecalibration(iLumi,cutoff,hb_recalib);
+    }
     if(hf_recalib)  hf_recalibration = new HFRecalibration();
 
     //     std::cout << " HcalHardcodeCalibrations:  iLumi = " <<  iLumi << std::endl;
@@ -359,16 +364,22 @@ std::auto_ptr<HcalRespCorrs> HcalHardcodeCalibrations::produceRespCorrs (const H
   rcd.getRecord<HcalRecNumberingRecord>().get(htopo);
   const HcalTopology* topo=&(*htopo);
 
-  //set depth segmentation for HE recalib - only happens once
-  if(he_recalibration && !setHEdsegm){
+  //set depth segmentation for HB/HE recalib - only happens once
+  if((he_recalibration && !setHEdsegm) || (hb_recalibration && !setHBdsegm)){
     std::vector<std::vector<int>> m_segmentation;
-    m_segmentation.resize(29);
-    for (int i = 0; i < 29; i++) {
+    m_segmentation.resize(30);
+    for (int i = 0; i < 30; i++) {
       if(i>0) topo->getDepthSegmentation(i,m_segmentation[i]);
     }
 
-    he_recalibration->setDsegm(m_segmentation);
-	setHEdsegm = true;
+	if(he_recalibration && !setHEdsegm){
+      he_recalibration->setDsegm(m_segmentation);
+	  setHEdsegm = true;
+	}
+	if(hb_recalibration && !setHBdsegm){
+      hb_recalibration->setDsegm(m_segmentation);
+	  setHBdsegm = true;
+	}
   }
 
   std::auto_ptr<HcalRespCorrs> result (new HcalRespCorrs (topo));
@@ -389,6 +400,19 @@ std::auto_ptr<HcalRespCorrs> HcalHardcodeCalibrations::produceRespCorrs (const H
 	<< "   corr = "  << corr << std::endl;
       */
 
+    }
+	else if ((hb_recalibration != 0 ) && 
+	((*cell).genericSubdet() == HcalGenericDetId::HcalGenBarrel)) {
+
+      int depth_ = HcalDetId(*cell).depth();
+      int ieta_  = HcalDetId(*cell).ieta();
+      corr = hb_recalibration->getCorr(ieta_, depth_); 
+
+      /*
+	std::cout << "HB ieta, depth = " << ieta_  << ",  " << depth_  
+	<< "   corr = "  << corr << std::endl;
+      */
+
     }
     else if ((hf_recalibration != 0 ) && 
 	((*cell).genericSubdet() == HcalGenericDetId::HcalGenForward)) {   

CalibCalorimetry/HcalPlugins/src/HcalHardcodeCalibrations.h

@@ -13,6 +13,7 @@
 #include "Geometry/CaloTopology/interface/HcalTopology.h"
 #include "CondFormats/HcalObjects/interface/AllObjects.h"
 #include "HERecalibration.h"
+#include "HBRecalibration.h"
 #include "DataFormats/HcalCalibObjects/interface/HFRecalibration.h"
 
 class ParameterSet;
@@ -89,8 +90,10 @@ class HcalHardcodeCalibrations : public edm::ESProducer,
 private:
   double iLumi;
   HERecalibration* he_recalibration;  
+  HBRecalibration* hb_recalibration;
   HFRecalibration* hf_recalibration;  
   bool switchGainWidthsForTrigPrims; 
   bool setHEdsegm;
+  bool setHBdsegm;
 };
 

DataFormats/HcalCalibObjects/interface/HBDarkening.h

@@ -0,0 +1,36 @@
+#ifndef HcalCalibObjects_HBDarkening_h
+#define HcalCalibObjects_HBDarkening_h
+//
+// Simple class with parameterizaed function to get darkening attenuation 
+// coefficiant for SLHC conditions
+// intlumi is integrated luminosity (fb-1),
+// layer is HB layer number (from 0 up to 16), NB: 1-17 in HcalTestNumbering
+//
+// scenario = 0 - full replacement of HB scintillators
+// scenario = 1 - no replacement, full stage darkening
+//
+
+class HBDarkening {
+
+public:
+  HBDarkening(unsigned int scenario);
+  ~HBDarkening();
+
+  float degradation(float intlumi, int ieta, int lay) const;
+
+  static const char* scenarioDescription (unsigned int scenario);
+
+  // 2 contsants below are used in CalibCalorimetry/HcalPlugins HERecalibration
+  // (1) number of HE ieta bins for darkening   
+  static const unsigned int nEtaBins = 16; 
+  // (2) max. number of HE scint. layers
+  static const unsigned int nScintLayers = 17;
+
+private:
+  int ieta_shift;
+  float lumiscale[nEtaBins][nScintLayers];
+
+};
+
+
+#endif // HBDarkening_h

DataFormats/HcalCalibObjects/src/HBDarkening.cc

@@ -0,0 +1,81 @@
+///////////////////////////////////////////////////////////////////////////////
+// File: HBDarkening.cc
+// Description: simple helper class containing parameterized function 
+//              to be used for the SLHC darkening calculation in HB 
+///////////////////////////////////////////////////////////////////////////////
+
+#include "DataFormats/HcalCalibObjects/interface/HBDarkening.h"
+#include <cmath>
+#include "FWCore/MessageLogger/interface/MessageLogger.h"
+
+//#define DebugLog
+
+HBDarkening::HBDarkening(unsigned int scenario) {
+  //HB starts at tower 1
+  ieta_shift = 1;
+
+  //scale parameter in fb-1 for exponential darkening from radiation damage for each tile
+  //extrapolated from HE radiation damage model (HEDarkening) using FLUKA dose map
+  //notes: extrapolation ignores geometrical and material differences in HB vs. HE
+  //       which affects the material averaging in the coarse FLUKA dose map
+  //       also, atmosphere in the barrel has more N2 and less O2 vs. HE
+  //conclusion: very rough extrapolation, caveat emptor!
+  float _lumiscale[nEtaBins][nScintLayers] = {
+      {887.1,887.1,1199.3,1314.0,1615.8,1954.0,2504.1,2561.7,3249.6,4553.7,5626.1,5725.6,6777.4,8269.8,10061.8,15253.3,22200.9},
+      {990.8,990.8,1140.8,1325.4,1676.9,2036.6,2900.5,2789.1,3460.2,4405.9,4184.0,5794.6,6157.4,7646.0,11116.1,18413.9,26813.3},
+      {971.2,971.2,1244.3,1456.2,1760.1,2299.1,2603.2,3012.3,3933.9,4787.3,4503.9,6624.3,7059.4,9369.5,12038.0,20048.8,23541.3},
+      {877.3,877.3,1145.2,1322.5,1604.9,1924.0,2893.6,2827.4,4085.0,5320.2,4740.6,5693.5,5715.1,7373.4,8305.1,16079.9,21702.3},
+      {919.8,919.8,1223.2,1376.9,1742.6,1964.7,2494.7,3335.6,4520.6,4869.5,4895.6,5740.4,7089.8,8765.9,10045.7,17408.0,24726.7},
+      {901.1,901.1,1114.3,1391.2,1733.2,2210.7,2733.8,3399.4,3715.2,3626.3,3371.3,4653.8,5911.6,7204.2,7584.7,4760.1,11156.8},
+      {904.0,904.0,1112.9,1351.9,1722.3,2008.2,2709.8,3101.9,3470.5,4679.0,5843.6,6343.8,7883.3,11266.8,16607.2,10882.3,25428.4},
+      {930.7,930.7,1225.3,1341.9,1744.0,2253.8,2805.1,3329.9,3665.6,5179.6,5677.8,5753.0,5662.3,9516.5,10769.4,13892.9,16661.1},
+      {953.2,953.2,1240.4,1487.3,1719.8,1930.6,2595.7,3172.5,3881.0,5247.5,4934.0,6576.4,6353.7,9259.2,12264.5,13261.8,12222.6},
+      {877.4,877.4,1114.1,1346.0,1604.9,1997.6,2708.9,3247.9,3704.4,4568.2,4984.4,7000.8,7896.7,7970.0,12555.6,10062.1,18386.8},
+      {876.2,876.2,1127.1,1336.0,1753.1,1944.4,2641.1,3445.1,3810.2,4033.0,5301.2,5170.7,6062.3,9815.8,12854.2,14297.3,20692.1},
+      {841.2,841.2,1051.1,1229.5,1576.5,1983.2,2282.2,2981.2,3271.7,4417.1,3765.2,4491.8,4626.6,7173.2,12953.0,7861.2,19338.6},
+      {848.2,848.2,1072.8,1228.5,1497.9,1876.1,2279.7,2744.2,3325.9,4021.8,4081.1,3750.6,4584.2,6170.8,9020.5,12058.8,13492.2},
+      {780.3,780.3,977.6,1103.2,1323.2,1548.3,1970.1,2217.5,2761.0,3049.0,2913.0,3832.3,4268.6,5242.1,6432.8,5999.5,6973.9},
+      {515.2,515.2,662.1,836.5,1110.4,1214.4,1664.4,1919.0,2341.0,2405.2,2647.5,2593.8,2586.3,2814.4,2826.8,0.0,0.0},
+      {409.3,409.3,489.2,700.3,960.2,1103.5,909.8,934.6,1148.6,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0}
+  };
+
+  //store array
+  //no factor needed to account for increased flux at 14 TeV vs 7-8 TeV
+  //already included in extrapolation from HE
+  for(unsigned int j = 0; j < nEtaBins; j++){
+    for(unsigned int i = 0; i < nScintLayers; i++){
+	  if (scenario == 0) lumiscale[j][i] = 0;
+      else lumiscale[j][i] = _lumiscale[j][i];
+	}
+  } 
+}
+
+HBDarkening::~HBDarkening() { }
+
+float HBDarkening::degradation(float intlumi, int ieta, int lay) const {
+  //no lumi, no darkening
+  if(intlumi <= 0) return 1.;
+
+  //shift ieta tower index
+  ieta -= ieta_shift;
+
+  //if outside eta range, no darkening
+  if(ieta < 0 || ieta >= (int)nEtaBins) return 1.;
+
+  //layer index does not need a shift in HB
+
+  //if outside layer range, no darkening
+  if(lay < 0 || lay >= (int)nScintLayers) return 1.;
+
+  // if replaced tile by scenario
+  if (lumiscale[ieta][lay] == 0.0) return 1.;
+
+  //return darkening factor
+  return (exp(-intlumi/lumiscale[ieta][lay]));
+}
+
+const char* HBDarkening::scenarioDescription (unsigned int scenario) {
+  if (scenario == 0) return "full replacement of HE scintillators, no darkening";
+  else if (scenario == 1) return "no replacement, full stage darkening";
+  return "undefined scenario: assume no replacement, full stage darkening";
+}