EnergyResolutionVsLumi.cc

Go to the documentation of this file.

#include <string>
#include <vector>
#include "SimG4CMS/Calo/interface/EvolutionECAL.h"
#include "SimG4CMS/Calo/interface/EnergyResolutionVsLumi.h"
#include "DataFormats/EcalDetId/interface/EBDetId.h"
#include "DataFormats/EcalDetId/interface/EEDetId.h"
#include "DataFormats/EcalDetId/interface/EcalSubdetector.h"

EnergyResolutionVsLumi::EnergyResolutionVsLumi()
{
  m_lumi=0;
  m_instlumi=0;  
}

EnergyResolutionVsLumi::~EnergyResolutionVsLumi()
{
}

EnergyResolutionVsLumi::DegradationAtEta EnergyResolutionVsLumi::CalculateDegradation(double eta)
{
  DegradationAtEta result;

  result.eta = eta;
  double totLumi = m_lumi;
  double instLumi = m_instlumi;

  EvolutionECAL model;

  // Index of induced absorption due to EM damages in PWO4
  result.muEM = model.InducedAbsorptionEM(instLumi, eta);
  
  // Index of induced absorption due to hadron damages in PWO4
  result.muHD = model.InducedAbsorptionHadronic(totLumi, eta);

  // Average degradation of amplitude due to transparency change
  result.ampDropTransparency = model.DegradationMeanEM50GeV(result.muEM+result.muHD);

  // Average degradation of amplitude due to photo-detector aging
  result.ampDropPhotoDetector = model.AgingVPT(instLumi, totLumi, eta);

  result.ampDropTotal = result.ampDropTransparency*result.ampDropPhotoDetector;

  // Noise increase in ADC counts due to photo-detector and front-end
  result.noiseIncreaseADC =  model.NoiseFactorFE(totLumi, eta);

  // Resolution degradation due to LY non-uniformity caused by transparency loss
  result.resolutitonConstantTerm =  model.ResolutionConstantTermEM50GeV(result.muEM+result.muHD);

  return result;
}

double EnergyResolutionVsLumi::calcmuEM(double eta)
{
  double instLumi = m_instlumi;
  EvolutionECAL model;
  double result = model.InducedAbsorptionEM(instLumi, eta);
  return result;
}

double EnergyResolutionVsLumi::calcmuHD(double eta)
{
  double totLumi = m_lumi;
  EvolutionECAL model;
  double result = model.InducedAbsorptionHadronic(totLumi, eta);
  return result;
}

void  EnergyResolutionVsLumi::calcmuTot(){

  for(int iEta=1; iEta<=EBDetId::MAX_IETA ;++iEta) {
    if(iEta==0) continue;

      double eta=EBDetId::approxEta(EBDetId(iEta,1));
      eta = std::abs(eta);
      double r= calcmuTot(eta);
      
      mu_eta[iEta]=r;
      vpt_eta[iEta]=1.0;

  }

  for(int iX=EEDetId::IX_MIN; iX<=EEDetId::IX_MAX ;++iX) {
    for(int iY=EEDetId::IY_MIN; iY<=EEDetId::IY_MAX; ++iY) {
      if (EEDetId::validDetId(iX,iY,1))
    {
      EEDetId eedetidpos(iX,iY,1);
      double eta= -log(tan(0.5*atan(sqrt((iX-50.0)*(iX-50.0)+(iY-50.0)*(iY-50.0))*2.98/328.)));
          eta = std::abs(eta);
          double r=calcmuTot(eta);
          double v=calcampDropPhotoDetector(eta);

      mu_eta[EBDetId::MAX_IETA+iX+iY*(EEDetId::IX_MAX)]=r;
      vpt_eta[EBDetId::MAX_IETA+iX+iY*(EEDetId::IX_MAX)]=v;
      //std::cout<<"eta/mu/vpt"<<eta<<"/"<<r<<"/"<<v<<std::endl;
    }
    }
  }
}

double EnergyResolutionVsLumi::calcLightCollectionEfficiencyWeighted(DetId id, double z)
{

  double v=1.0;
  double muTot=0;
  if(id.subdetId()==EcalBarrel) {
    EBDetId ebId(id);
    int ieta= std::abs(ebId.ieta());
    muTot= mu_eta[ieta];
    
  } else if(id.subdetId()==EcalEndcap){
    EEDetId eeId(id);
    int ix= eeId.ix();
    int iy= eeId.iy();
    
    muTot= mu_eta[EBDetId::MAX_IETA+ix+iy*(EEDetId::IX_MAX)];
    v=vpt_eta[EBDetId::MAX_IETA+ix+iy*(EEDetId::IX_MAX)];
  } else {
    muTot=0;
  }
  double zcor=z;
  EvolutionECAL model;
  if(z<0.02 )     { zcor=0.02; }
  else if(z>0.98) { zcor=0.98; }

  double result=model.LightCollectionEfficiencyWeighted( zcor , muTot)*v;
  
  return result; 
}

double EnergyResolutionVsLumi::calcLightCollectionEfficiencyWeighted2(double eta, double z, double mu_ind)
{
  if(mu_ind<0) mu_ind=this->calcmuTot(eta);
  double v= this->calcampDropPhotoDetector(eta);
  EvolutionECAL model;
  double result=model.LightCollectionEfficiencyWeighted( z , mu_ind)*v;
  return result; 
}

double EnergyResolutionVsLumi::calcmuTot(double eta)
{
  double totLumi = m_lumi;
  double instLumi = m_instlumi;
  EvolutionECAL model;
  double muEM = model.InducedAbsorptionEM(instLumi, eta);
  double muH = model.InducedAbsorptionHadronic(totLumi, eta);
  double result=muEM+muH;
  return result;
} 

double EnergyResolutionVsLumi::calcampDropTransparency(double eta)
{
  double muEM=this->calcmuEM(eta);
  double muHD=this->calcmuHD(eta);
  EvolutionECAL model;
  double result = model.DegradationMeanEM50GeV(muEM+muHD);
  return result;
}

double EnergyResolutionVsLumi::calcampDropPhotoDetector(double eta)
{
  double instLumi = m_instlumi;
  double totLumi = m_lumi;
  EvolutionECAL model;
  double result = model.AgingVPT(instLumi, totLumi, eta);
  return result;
}

double EnergyResolutionVsLumi::calcampDropTotal(double eta)
{
  double tra= this->calcampDropTransparency(eta);
  double pho= this->calcampDropPhotoDetector(eta);
  double result = tra*pho;
  return result;
}

double EnergyResolutionVsLumi::calcnoiseIncreaseADC(double eta)
{
  double totLumi = m_lumi;
  EvolutionECAL model;
  double result = model.NoiseFactorFE(totLumi, eta);
  return result;
  // noise increase in ADC
}

double EnergyResolutionVsLumi::calcnoiseADC(double eta)
{
  double totLumi = m_lumi;
  double Nadc = 1.1;
  double result =1.0;
  EvolutionECAL model;
  if(std::abs(eta)<1.497){
    Nadc = 1.1;
    result = model.NoiseFactorFE(totLumi, eta)*Nadc;
  }else{
    Nadc = 2.0;
    result=Nadc;
    // endcaps no increase in ADC
  }
  return result;

}

double EnergyResolutionVsLumi::calcresolutitonConstantTerm(double eta)
{
  double muEM=this->calcmuEM(eta);
  double muHD=this->calcmuHD(eta);
  EvolutionECAL model;
  double result = model.ResolutionConstantTermEM50GeV(muEM+muHD);
  return result;
}

double  EnergyResolutionVsLumi::Resolution(double eta, double ene)
{  
  // Initial parameters for ECAL resolution
  double S;
  double Nadc;
  double adc2GeV;
  double C;
  if(eta<1.497){
    S = 0.028;           // CMS note 2006/148 (EB test beam)
    Nadc = 1.1; 
    adc2GeV = 0.039;
    C = 0.003;
  }else{
    S = 0.052;          //  CMS DN 2009/002
    Nadc = 2.0;
    adc2GeV = 0.069;
    C = 0.004;
  }

  DegradationAtEta d = CalculateDegradation(eta);

  // adjust resolution parameters
  S /= sqrt(d.ampDropTotal);
  Nadc *= d.noiseIncreaseADC;
  adc2GeV /= d.ampDropTotal;
  double N = Nadc*adc2GeV*3.0;   // 3x3 noise in GeV 
  C = sqrt(C*C + d.resolutitonConstantTerm*d.resolutitonConstantTerm);

  return sqrt(S*S/ene + N*N/ene/ene + C*C);

}
/*
void EnergyResolutionVsLumi::Decomposition()
{
  double eta = 2.2;
  m_instlumi = 5.0e+34;
  m_lumi = 3000.0;

  DegradationAtEta d = this->CalculateDegradation(eta);

  // Initial parameters for ECAL resolution
  double S;
  double Nadc;
  double adc2GeV;
  double C;
  if(eta<1.497){
    S = 0.028;           // CMS note 2006/148 (EB test beam)
    Nadc = 1.1; 
    adc2GeV = 0.039;
    C = 0.003;
  }else{
    S = 0.052;          //  CMS DN 2009/002
    Nadc = 2.0;
    adc2GeV = 0.069;
    C = 0.0038;
  }

  // adjust resolution parameters
  S /= sqrt(d.ampDropTotal);
  Nadc *= d.noiseIncreaseADC;
  adc2GeV /= d.ampDropTotal;
  //  double N = Nadc*adc2GeV*3.0;   // 3x3 noise in GeV 
  C = sqrt(C*C + d.resolutitonConstantTerm*d.resolutitonConstantTerm);

  for(double ene=1.0; ene<1e+3; ene*=1.1){
    // this is the resolution
    double res =  sqrt(S*S/ene + N*N/ene/ene + C*C);
    double factor = 1.0;
    factor = sin(2.0*atan(exp(-1.0*eta)));
  }
}
*/