Calo/src/EnergyResolutionVsLumi.cc

0001 #include <string>
0002 #include <vector>
0003 #include "SimG4CMS/Calo/interface/EvolutionECAL.h"
0004 #include "SimG4CMS/Calo/interface/EnergyResolutionVsLumi.h"
0005 #include "DataFormats/EcalDetId/interface/EBDetId.h"
0006 #include "DataFormats/EcalDetId/interface/EEDetId.h"
0007 #include "DataFormats/EcalDetId/interface/EcalSubdetector.h"
0008
0009 EnergyResolutionVsLumi::EnergyResolutionVsLumi() {
0010   m_lumi = 0;
0011   m_instlumi = 0;
0012 }
0013
0014 EnergyResolutionVsLumi::DegradationAtEta EnergyResolutionVsLumi::CalculateDegradation(double eta) {
0015   DegradationAtEta result;
0016
0017   result.eta = eta;
0018   double totLumi = m_lumi;
0019   double instLumi = m_instlumi;
0020
0021   EvolutionECAL model;
0022
0023   // Index of induced absorption due to EM damages in PWO4
0024   result.muEM = model.InducedAbsorptionEM(instLumi, eta);
0025
0026   // Index of induced absorption due to hadron damages in PWO4
0027   result.muHD = model.InducedAbsorptionHadronic(totLumi, eta);
0028
0029   // Average degradation of amplitude due to transparency change
0030   result.ampDropTransparency = model.DegradationMeanEM50GeV(result.muEM + result.muHD);
0031
0032   // Average degradation of amplitude due to photo-detector aging
0033   result.ampDropPhotoDetector = model.AgingVPT(instLumi, totLumi, eta);
0034
0035   result.ampDropTotal = result.ampDropTransparency * result.ampDropPhotoDetector;
0036
0037   // Noise increase in ADC counts due to photo-detector and front-end
0038   result.noiseIncreaseADC = model.NoiseFactorFE(totLumi, eta);
0039
0040   // Resolution degradation due to LY non-uniformity caused by transparency loss
0041   result.resolutitonConstantTerm = model.ResolutionConstantTermEM50GeV(result.muEM + result.muHD);
0042
0043   return result;
0044 }
0045
0046 double EnergyResolutionVsLumi::calcmuEM(double eta) {
0047   double instLumi = m_instlumi;
0048   EvolutionECAL model;
0049   double result = model.InducedAbsorptionEM(instLumi, eta);
0050   return result;
0051 }
0052
0053 double EnergyResolutionVsLumi::calcmuHD(double eta) {
0054   double totLumi = m_lumi;
0055   EvolutionECAL model;
0056   double result = model.InducedAbsorptionHadronic(totLumi, eta);
0057   return result;
0058 }
0059
0060 void EnergyResolutionVsLumi::calcmuTot() {
0061   for (int iEta = 1; iEta <= EBDetId::MAX_IETA; ++iEta) {
0062     if (iEta == 0)
0063       continue;
0064
0065     double eta = EBDetId::approxEta(EBDetId(iEta, 1));
0066     eta = std::abs(eta);
0067     double r = calcmuTot(eta);
0068
0069     mu_eta[iEta] = r;
0070     vpt_eta[iEta] = 1.0;
0071   }
0072
0073   for (int iX = EEDetId::IX_MIN; iX <= EEDetId::IX_MAX; ++iX) {
0074     for (int iY = EEDetId::IY_MIN; iY <= EEDetId::IY_MAX; ++iY) {
0075       if (EEDetId::validDetId(iX, iY, 1)) {
0076         EEDetId eedetidpos(iX, iY, 1);
0077         double eta = -log(tan(0.5 * atan(sqrt((iX - 50.0) * (iX - 50.0) + (iY - 50.0) * (iY - 50.0)) * 2.98 / 328.)));
0078         eta = std::abs(eta);
0079         double r = calcmuTot(eta);
0080         double v = calcampDropPhotoDetector(eta);
0081
0082         mu_eta[EBDetId::MAX_IETA + iX + iY * (EEDetId::IX_MAX)] = r;
0083         vpt_eta[EBDetId::MAX_IETA + iX + iY * (EEDetId::IX_MAX)] = v;
0084         //std::cout<<"eta/mu/vpt"<<eta<<"/"<<r<<"/"<<v<<std::endl;
0085       }
0086     }
0087   }
0088 }
0089
0090 double EnergyResolutionVsLumi::calcLightCollectionEfficiencyWeighted(DetId id, double z) {
0091   double v = 1.0;
0092   double muTot = 0;
0093   if (id.subdetId() == EcalBarrel) {
0094     EBDetId ebId(id);
0095     int ieta = std::abs(ebId.ieta());
0096     muTot = mu_eta[ieta];
0097
0098   } else if (id.subdetId() == EcalEndcap) {
0099     EEDetId eeId(id);
0100     int ix = eeId.ix();
0101     int iy = eeId.iy();
0102
0103     muTot = mu_eta[EBDetId::MAX_IETA + ix + iy * (EEDetId::IX_MAX)];
0104     v = vpt_eta[EBDetId::MAX_IETA + ix + iy * (EEDetId::IX_MAX)];
0105   } else {
0106     muTot = 0;
0107   }
0108   double zcor = z;
0109   EvolutionECAL model;
0110   if (z < 0.02) {
0111     zcor = 0.02;
0112   } else if (z > 0.98) {
0113     zcor = 0.98;
0114   }
0115
0116   double result = model.LightCollectionEfficiencyWeighted(zcor, muTot) * v;
0117
0118   return result;
0119 }
0120
0121 double EnergyResolutionVsLumi::calcLightCollectionEfficiencyWeighted2(double eta, double z, double mu_ind) {
0122   if (mu_ind < 0)
0123     mu_ind = this->calcmuTot(eta);
0124   double v = this->calcampDropPhotoDetector(eta);
0125   EvolutionECAL model;
0126   double result = model.LightCollectionEfficiencyWeighted(z, mu_ind) * v;
0127   return result;
0128 }
0129
0130 double EnergyResolutionVsLumi::calcmuTot(double eta) {
0131   double totLumi = m_lumi;
0132   double instLumi = m_instlumi;
0133   EvolutionECAL model;
0134   double muEM = model.InducedAbsorptionEM(instLumi, eta);
0135   double muH = model.InducedAbsorptionHadronic(totLumi, eta);
0136   double result = muEM + muH;
0137   return result;
0138 }
0139
0140 double EnergyResolutionVsLumi::calcampDropTransparency(double eta) {
0141   double muEM = this->calcmuEM(eta);
0142   double muHD = this->calcmuHD(eta);
0143   EvolutionECAL model;
0144   double result = model.DegradationMeanEM50GeV(muEM + muHD);
0145   return result;
0146 }
0147
0148 double EnergyResolutionVsLumi::calcampDropPhotoDetector(double eta) {
0149   double instLumi = m_instlumi;
0150   double totLumi = m_lumi;
0151   EvolutionECAL model;
0152   double result = model.AgingVPT(instLumi, totLumi, eta);
0153   return result;
0154 }
0155
0156 double EnergyResolutionVsLumi::calcampDropTotal(double eta) {
0157   double tra = this->calcampDropTransparency(eta);
0158   double pho = this->calcampDropPhotoDetector(eta);
0159   double result = tra * pho;
0160   return result;
0161 }
0162
0163 double EnergyResolutionVsLumi::calcnoiseIncreaseADC(double eta) {
0164   double totLumi = m_lumi;
0165   EvolutionECAL model;
0166   double result = model.NoiseFactorFE(totLumi, eta);
0167   return result;
0168   // noise increase in ADC
0169 }
0170
0171 double EnergyResolutionVsLumi::calcnoiseADC(double eta) {
0172   double totLumi = m_lumi;
0173   double Nadc = 1.1;
0174   double result = 1.0;
0175   EvolutionECAL model;
0176   if (std::abs(eta) < 1.497) {
0177     Nadc = 1.1;
0178     result = model.NoiseFactorFE(totLumi, eta) * Nadc;
0179   } else {
0180     Nadc = 2.0;
0181     result = Nadc;
0182     // endcaps no increase in ADC
0183   }
0184   return result;
0185 }
0186
0187 double EnergyResolutionVsLumi::calcresolutitonConstantTerm(double eta) {
0188   double muEM = this->calcmuEM(eta);
0189   double muHD = this->calcmuHD(eta);
0190   EvolutionECAL model;
0191   double result = model.ResolutionConstantTermEM50GeV(muEM + muHD);
0192   return result;
0193 }
0194
0195 double EnergyResolutionVsLumi::Resolution(double eta, double ene) {
0196   // Initial parameters for ECAL resolution
0197   double S;
0198   double Nadc;
0199   double adc2GeV;
0200   double C;
0201   if (eta < 1.497) {
0202     S = 0.028;  // CMS note 2006/148 (EB test beam)
0203     Nadc = 1.1;
0204     adc2GeV = 0.039;
0205     C = 0.003;
0206   } else {
0207     S = 0.052;  //  CMS DN 2009/002
0208     Nadc = 2.0;
0209     adc2GeV = 0.069;
0210     C = 0.004;
0211   }
0212
0213   DegradationAtEta d = CalculateDegradation(eta);
0214
0215   // adjust resolution parameters
0216   S /= sqrt(d.ampDropTotal);
0217   Nadc *= d.noiseIncreaseADC;
0218   adc2GeV /= d.ampDropTotal;
0219   double N = Nadc * adc2GeV * 3.0;  // 3x3 noise in GeV
0220   C = sqrt(C * C + d.resolutitonConstantTerm * d.resolutitonConstantTerm);
0221
0222   return sqrt(S * S / ene + N * N / ene / ene + C * C);
0223 }
0224 /*
0225 void EnergyResolutionVsLumi::Decomposition()
0226 {
0227   double eta = 2.2;
0228   m_instlumi = 5.0e+34;
0229   m_lumi = 3000.0;
0230
0231   DegradationAtEta d = this->CalculateDegradation(eta);
0232
0233   // Initial parameters for ECAL resolution
0234   double S;
0235   double Nadc;
0236   double adc2GeV;
0237   double C;
0238   if(eta<1.497){
0239     S = 0.028;           // CMS note 2006/148 (EB test beam)
0240     Nadc = 1.1;
0241     adc2GeV = 0.039;
0242     C = 0.003;
0243   }else{
0244     S = 0.052;          //  CMS DN 2009/002
0245     Nadc = 2.0;
0246     adc2GeV = 0.069;
0247     C = 0.0038;
0248   }
0249
0250   // adjust resolution parameters
0251   S /= sqrt(d.ampDropTotal);
0252   Nadc *= d.noiseIncreaseADC;
0253   adc2GeV /= d.ampDropTotal;
0254   //  double N = Nadc*adc2GeV*3.0;   // 3x3 noise in GeV
0255   C = sqrt(C*C + d.resolutitonConstantTerm*d.resolutitonConstantTerm);
0256
0257   for(double ene=1.0; ene<1e+3; ene*=1.1){
0258     // this is the resolution
0259     double res =  sqrt(S*S/ene + N*N/ene/ene + C*C);
0260     double factor = 1.0;
0261     factor = sin(2.0*atan(exp(-1.0*eta)));
0262   }
0263 }
0264 */