Skip to content

Commit

Permalink
latest calibrations from Pedro + global compensation plugin
Browse files Browse the repository at this point in the history
  • Loading branch information
@lgray
lgray committed Mar 15, 2015
1 parent e739ebf commit 3b215b1
Show file tree
Hide file tree
Showing 4 changed files with 224 additions and 49 deletions.
43 changes: 42 additions & 1 deletion RecoParticleFlow/PandoraTranslator/data/PandoraSettingsBasic_cms.xml
View file
Expand Up @@ -7,7 +7,48 @@
<ShouldCollapseMCParticlesToPfoTarget>true</ShouldCollapseMCParticlesToPfoTarget>

<!-- PLUGIN SETTINGS -->
<HadronicEnergyCorrectionPlugins>CleanClusters ScaleHotHadrons</HadronicEnergyCorrectionPlugins>
<HadronicEnergyCorrectionPlugins>CleanClusters CMSGlobalHadronCompensation ScaleHotHadrons</HadronicEnergyCorrectionPlugins>

<CMSGlobalHadronCompensation>
<MipEnergyThreshold>10.0</MipEnergyThreshold>
<IntegMIP_emEn_EE>0.423</IntegMIP_emEn_EE>
<!-- e_em_BH function: [0]*x -->
<e_em_BH> 4.14730994206</e_em_BH>
<!-- e_em_EE function: [0]*x -->
<e_em_EE> 4.41625022557 </e_em_EE>
<!-- e_em_FH function: [0]*x -->
<e_em_FH> 5.6226096912 </e_em_FH>

<!-- BY MEDIAN -->
<!-- pioe_BH function: ([0]*(((1-exp(-([1]*x)))/(1+exp(-([2]*x))))*exp(-([3]*x))))+[4] -->
<!-- <pioe_BH> 0.715561124085 0.201541211354 0.00679202388913 0.00163121938632 0.54252374821 </pioe_BH> -->
<!-- pioe_FH function: ([0]*(((1-exp(-([1]*x)))/(1+exp(-([2]*x))))*exp(-([3]*x))))+[4] -->
<!-- <pioe_FH> 0.430841020284 0.252699052283 0.0380850086469 0.0 0.50349718687 </pioe_FH> -->
<!-- pioe_EE function: ([0]*(((1-exp(-([1]*x)))/(1+exp(-([2]*x))))*exp(-([3]*x))))+[4] -->
<!-- <pioe_EE> 0.615988413399 9.99248772533 2.22098017005 0.0 0.380425357427 </pioe_EE> -->
<!-- pioe_bananaParam_1 function: [0]+([1]*(1-exp(-([2]*x)))) -->
<!-- <pioe_bananaParam_1> 0.165517562394 -0.489815099048 0.0388102416692 </pioe_bananaParam_1> -->
<!-- pioe_bananaParam_0 function: [0]+([1]*(1-exp(-([2]*x)))) -->
<!-- <pioe_bananaParam_0> 0.977178682989 0.0252398439787 0.26978497941 </pioe_bananaParam_0> -->
<!-- pioe_bananaParam_2 function: [0]+([1]*(1-exp(-([2]*x)))) -->
<!-- <pioe_bananaParam_2> -0.145929892662 0.447931176619 0.0490124900322 </pioe_bananaParam_2> -->

<!-- BY MODE -->
<!-- pioe_HEB function: ([0]*(((1-exp(-([1]*x)))/(1+exp(-([2]*x))))*exp(-([3]*x))))+[4] -->
<pioe_HEB> 0.842152103125 0.123113884222 0.00499766525242 0.00155314804086 0.500339972452 </pioe_HEB>
<!-- pioe_HEF function: ([0]*(((1-exp(-([1]*x)))/(1+exp(-([2]*x))))*exp(-([3]*x))))+[4] -->
<pioe_HEF> 0.415875311091 0.1651386381 0.0351926420879 0.0 0.51838818602 </pioe_HEF>
<!-- pioe_EE function: ([0]*(((1-exp(-([1]*x)))/(1+exp(-([2]*x))))*exp(-([3]*x))))+[4] -->
<pioe_EE> 0.58304638498 9.08251124371 2.21320527452 0.0 0.413017771562 </pioe_EE>
<!-- pioe_bananaParam_1 function: [0]+([1]*(1-exp(-([2]*x)))) -->
<pioe_bananaParam_1> 0.192256003375 -0.515170751706 0.0384348463059 </pioe_bananaParam_1>
<!-- pioe_bananaParam_0 function: [0]+([1]*(1-exp(-([2]*x)))) -->
<pioe_bananaParam_0> 0.975281059636 0.0266241561011 0.282004059773 </pioe_bananaParam_0>
<!-- pioe_bananaParam_2 function: [0]+([1]*(1-exp(-([2]*x)))) -->
<pioe_bananaParam_2> -0.109883036753 0.409838061737 0.0465469104675 </pioe_bananaParam_2>

</CMSGlobalHadronCompensation>

<EmShowerPlugin>
LCEmShowerId
<RmsCutEnergy>70.0</RmsCutEnergy>
Expand Down
54 changes: 54 additions & 0 deletions RecoParticleFlow/PandoraTranslator/data/pandoraCalibrPars_pedro05032015.txt
View file
@@ -0,0 +1,54 @@
Calibr_ADC2GeV_EE 1.
Calibr_ADC2GeV_HEF 1.
Calibr_ADC2GeV_HEB 1.

#additional calibration
#this number is energy scale/MIP from Pedro's slides
EMaddCalibrEE 219.60072595
EMaddCalibrHEF 191.76470588
EMaddCalibrHEB 14.949279231

#mu_k * alpha_k (calibrated to lambdas) /MIP from Pedro's slides
# slope * mu_i * alpha_i_lambda / MIP-in-GeV
HADaddCalibrEE 4676.8435572
HADaddCalibrHEF 2360.6841204
HADaddCalibrHEB 176.1126

# #threshold in GeV -- increase a little bit
ECalThresEndCap 27.75e-6
HCalThresEndCapHEF 42.70e-6
HCalThresEndCapHEB 760.2e-6
ECalThresBarrel 0.
HCalThresBarrel 0.

ECalMipThresEndCap 0.5
HCalMipThresEndCapHEF 0.5
HCalMipThresEndCapHEB 0.5
ECalMipThresBarrel 10.
HCalMipThresBarrel 10.

# mip calibration ----------- 1/MIP-in-GeV
ECalToMipEndCap 18148.82
HCalToMipEndCapHEF 11764.71
HCalToMipEndCapHEB 667.38
HCalToMipBarrel 3.333333
ECalToMipBarrel 3.333333

MuonToMip 10.0

# EM calibration ----------------
ECalToEMGeVEndCap 1.0
HCalToEMGeVEndCapHEF 1.0
HCalToEMGeVEndCapHEB 1.0
ECalToEMGeVBarrel 1.0
HCalToEMGeVBarrel 1.0

# Hadron calibration
ECalToHadGeVEndCap 1.0
HCalToHadGeVEndCapHEF 1.0
HCalToHadGeVEndCapHEB 1.0
ECalToHadGeVBarrel 1.0
HCalToHadGeVBarrel 1.0



2 changes: 1 addition & 1 deletion RecoParticleFlow/PandoraTranslator/python/runPandora_cfi.py
View file
Expand Up @@ -17,7 +17,7 @@
# energyCorrMethod = cms.string('WEIGHTING'),
energyWeightFile = cms.FileInPath('RecoParticleFlow/PandoraTranslator/data/energyWeight.txt'),

calibrParFile = cms.FileInPath('RecoParticleFlow/PandoraTranslator/data/pandoraCalibrPars_pedro24022015.txt'),
calibrParFile = cms.FileInPath('RecoParticleFlow/PandoraTranslator/data/pandoraCalibrPars_pedro05032015.txt'),
layerDepthFile = cms.FileInPath('RecoParticleFlow/PandoraTranslator/data/HGCmaterial_v5.root'),
overburdenDepthFile = cms.FileInPath('RecoParticleFlow/PFClusterProducer/data/HGCMaterialOverburden.root'),
useOverburdenCorrection = cms.bool(False), #disabled until the overburden values make sense
Expand Down
174 changes: 127 additions & 47 deletions RecoParticleFlow/PandoraTranslator/src/CMSGlobalHadronCompensationPlugin.cc
View file
Expand Up @@ -12,16 +12,18 @@

#include "RecoParticleFlow/PandoraTranslator/interface/CMSGlobalHadronCompensationPlugin.h"

#include <DataFormats/ForwardDetId/interface/ForwardSubdetector.h>
#include <DataFormats/ParticleFlowReco/interface/PFRecHit.h>

#include <algorithm>

using namespace pandora;
using namespace cms_content;

//------------------------------------------------------------------------------------------------------------------------------------------

GlobalHadronCompensation::GlobalHadronCompensation() :
m_nMIPsCut(10.0f)
m_nMIPsCut(1.0),
m_e_em_EE(1.0), m_e_em_FH(1.0), m_e_em_BH(1.0)
{
}

Expand All @@ -36,53 +38,111 @@ StatusCode GlobalHadronCompensation::MakeEnergyCorrections(const Cluster *const

OrderedCaloHitList::const_iterator layer = orderedCaloHitList.begin();
OrderedCaloHitList::const_iterator hits_end = orderedCaloHitList.end();

float en_avg(0.f), nHits(0.f);

for( ; layer != hits_end; ++ layer ) {
const CaloHitList& hits_in_layer = *(layer->second);
for( const auto& hit : hits_in_layer ) {
// hack so that we can know if we are in the HEF or not (go back to cmssw det id)
const void* void_hit_ptr = hit->GetParentCaloHitAddress();
const reco::PFRecHit* original_hit_ptr = static_cast<const reco::PFRecHit*>(void_hit_ptr);
const uint32_t rawid = original_hit_ptr->detId();
const int subDetId = (rawid>>25)&0x7;
if( subDetId != ForwardSubdetector::HGCHEF ) continue;

//compute raw energy
float rawen_EE(0.f), rawen_FH(0.f), rawen_BH(0.f);
float avgen_mip_FH(0.f), nHits_FH(0.f);
for( ; layer != hits_end; ++ layer )
{
const CaloHitList& hits_in_layer = *(layer->second);
for( const auto& hit : hits_in_layer )
{
//check sub-detector id
const void* void_hit_ptr = hit->GetParentCaloHitAddress();
const reco::PFRecHit* original_hit_ptr = static_cast<const reco::PFRecHit*>(void_hit_ptr);
const uint32_t rawid = original_hit_ptr->detId();
const int subDetId = (rawid>>25)&0x7;

const float en_mip = hit->GetMipEquivalentEnergy();
en_avg += en_mip;
nHits += 1.f;
}
}
float nLambdas(hit->GetNCellInteractionLengths());
float nMIPs(hit->GetMipEquivalentEnergy());

if(subDetId == ForwardSubdetector::HGCEE)
{
rawen_EE += nLambdas*nMIPs/m_e_em_EE;
}
else if(subDetId == ForwardSubdetector::HGCHEB)
{
rawen_BH += nLambdas*nMIPs/m_e_em_FH;
}
else if( subDetId == ForwardSubdetector::HGCHEF )
{
rawen_FH += nLambdas*nMIPs/m_e_em_BH;
avgen_mip_FH += nMIPs;
nHits_FH += 1.0;
}
}
}

if( nHits == 0.f ) {
correctedHadronicEnergy *= 1.f;
return STATUS_CODE_SUCCESS;
}

en_avg /= nHits;

float nHits_avg(0.f), nHits_elim(0.f);
for( ; layer != hits_end; ++ layer ) {
const CaloHitList& hits_in_layer = *(layer->second);
for( const auto& hit : hits_in_layer ) {
// hack so that we can know if we are in the HEF or not (go back to cmssw det id)
const void* void_hit_ptr = hit->GetParentCaloHitAddress();
const reco::PFRecHit* original_hit_ptr = static_cast<const reco::PFRecHit*>(void_hit_ptr);
const uint32_t rawid = original_hit_ptr->detId();
const int subDetId = (rawid>>25)&0x7;
if( subDetId != ForwardSubdetector::HGCHEF ) continue;

const float en_mip = hit->GetMipEquivalentEnergy();
if( en_mip > en_avg ) {
nHits_avg += 1.f;
} else {
nHits_elim += 1.f;
}
//further correction only makes sense if enough energy
float rawen(rawen_EE+rawen_FH+rawen_BH);
if( rawen<0.01f)
{
correctedHadronicEnergy *= 1.f;
return STATUS_CODE_SUCCESS;
}
}

correctedHadronicEnergy *= (nHits-nHits_elim)/(nHits-nHits_avg);

//global compensation if energy in FH is significant
float c_FH(1.0);
if(nHits_FH>0)
{
avgen_mip_FH/=nHits_FH;
layer = orderedCaloHitList.begin();
float nHits_avg_FH(0.0), nHits_elim_FH(0.0);
for( ; layer != hits_end; ++ layer )
{
const CaloHitList& hits_in_layer = *(layer->second);
for( const auto& hit : hits_in_layer )
{
// hack so that we can know if we are in the HEF or not (go back to cmssw det id)
const void* void_hit_ptr = hit->GetParentCaloHitAddress();
const reco::PFRecHit* original_hit_ptr = static_cast<const reco::PFRecHit*>(void_hit_ptr);
const uint32_t rawid = original_hit_ptr->detId();
const int subDetId = (rawid>>25)&0x7;
if( subDetId != ForwardSubdetector::HGCHEF ) continue;

const float nMIPs = hit->GetMipEquivalentEnergy();
nHits_avg_FH += ( nMIPs > avgen_mip_FH);
nHits_elim_FH += ( nMIPs > m_nMIPsCut );
}
}

c_FH=(nHits_FH-nHits_elim_FH)/(nHits_FH-nHits_avg_FH);
}


//pi/e corrected energy
float en_EE(getByPiOverECorrectedEn(rawen_EE,rawen,int(ForwardSubdetector::HGCEE)));
float en_FH(c_FH*getByPiOverECorrectedEn(rawen_FH,rawen,int(ForwardSubdetector::HGCHEF)));
float en_BH(getByPiOverECorrectedEn(rawen_BH,rawen,int(ForwardSubdetector::HGCHEB)));
float en_HCAL(en_FH+en_BH);
float en(en_EE+en_HCAL);

//energy sharing (ECAL MIP subtracted)
float en_m_EEmip(std::max(en_EE-m_IntegMIP_emEn_EE,0.f)+en_HCAL);
float enFracInHCAL_m_mip(-1);
if(en_HCAL==0) enFracInHCAL_m_mip=0;
if(en_m_EEmip>0) enFracInHCAL_m_mip=en_HCAL/en_m_EEmip;

//apply residual correction according to energy sharing
float residualScale=getResidualScale(en,enFracInHCAL_m_mip);


/*
std::cout << "Initial: " << correctedHadronicEnergy
<< " Raw: " << rawen
<< " After pi/e: " << en
<< " GC: " << c_FH
<< " Residual: " << residualScale
<< " Final: " << residualScale*en
<< std::endl
<< "\t EE: " << rawen_EE << "->" << en_EE
<< " FH: " << rawen_FH << "->" << en_FH
<< " BH: " << rawen_BH << "->" << en_BH
<< std::endl;
*/
//all done here
correctedHadronicEnergy=residualScale*en;

return STATUS_CODE_SUCCESS;
}
Expand Down Expand Up @@ -110,8 +170,28 @@ float GlobalHadronCompensation::GetHadronicEnergyInLayer(const OrderedCaloHitLis

StatusCode GlobalHadronCompensation::ReadSettings(const TiXmlHandle xmlHandle)
{
PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle,
"MipEnergyThreshold", m_nMIPsCut));
PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "MipEnergyThreshold", m_nMIPsCut));
PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "IntegMIP_emEn_EE", m_IntegMIP_emEn_EE));
PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "e_em_EE", m_e_em_EE));
PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "e_em_FH", m_e_em_FH));
PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "e_em_BH", m_e_em_BH));
PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadVectorOfValues(xmlHandle, "pioe_EE", m_pioe_EE));
PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadVectorOfValues(xmlHandle, "pioe_FH", m_pioe_FH));
PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadVectorOfValues(xmlHandle, "pioe_BH", m_pioe_BH));
PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadVectorOfValues(xmlHandle, "pioe_bananaParam_0", m_pioe_bananaParam_0 ));
PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadVectorOfValues(xmlHandle, "pioe_bananaParam_1", m_pioe_bananaParam_1 ));
PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadVectorOfValues(xmlHandle, "pioe_bananaParam_2", m_pioe_bananaParam_2 ));

m_IntegMIP_emEn_EE=getByPiOverECorrectedEn(m_IntegMIP_emEn_EE,m_IntegMIP_emEn_EE,(int)(ForwardSubdetector::HGCEE));

std::cout << "[CMSGlobalHadronCompensationPlugin::ReadSettings]" << std::endl
<< "MIP threshold for GC: " << m_nMIPsCut << std::endl
<< "MIP value in EE (after pi/e/):" << m_IntegMIP_emEn_EE << std::endl
<< "e.m. scales " << 1./m_e_em_EE << " " << 1./m_e_em_FH << " " << 1./m_e_em_BH << std::endl
<< "# pi/e parameters " << m_pioe_EE.size() << " " << m_pioe_FH.size() << " " << m_pioe_BH.size() <<" " << std::endl
<< "# res. corr parameters " << m_pioe_bananaParam_0.size() << " " << m_pioe_bananaParam_1.size()<< " " << m_pioe_bananaParam_2.size() << std::endl;



return STATUS_CODE_SUCCESS;
}

0 comments on commit 3b215b1

Please sign in to comment.