Skip to content
Snippets Groups Projects
Commit 152aec76 authored by Emma Torro Pastor's avatar Emma Torro Pastor Committed by Graeme Stewart
Browse files

Adding a fix to reduce CalRatio pileup dependency (ATR-15255) (TrigL2LongLivedParticles-00-00-58)

parent df0166ff
No related branches found
No related tags found
No related merge requests found
......@@ -32,6 +32,7 @@ class TrigJetSplitter : public HLT::AllTEAlgo {
double m_minJetEt;
double m_maxJetEta;
double m_logRatio;
double m_pufixLogRatio;
bool m_reversedCut;
......
......@@ -9,8 +9,14 @@ from AthenaCommon.GlobalFlags import globalflags
from AthenaCommon.AppMgr import ServiceMgr
def getJetSplitterInstance( ):
return JetSplitter( name="JetSplitter" )
def getJetSplitterInstance( instance, logratio, pufixlogratio ):
name=instance+"_"+str(int(pufixlogratio*10))
return JetSplitter( name=name,
logratio=logratio,
pufixlogratio=pufixlogratio
)
def getJetSplitterInstance_LowLogRatio( ):
return JetSplitter_LowLogRatio( name="JetSplitter_LowLogRatio" )
......@@ -19,15 +25,15 @@ def getJetSplitterInstance_LowLogRatio( ):
class JetSplitter (TrigJetSplitter):
__slots__ = []
def __init__(self, name):
def __init__(self, name, logratio, pufixlogratio):
super( JetSplitter, self ).__init__( name )
self.JetInputKey = "TrigJetRec"
self.JetOutputKey = "SplitJet"
self.EtaHalfWidth = 0.4
self.PhiHalfWidth = 0.4
self.JetLogRatio = 1.2
self.JetLogRatio = logratio
self.JetPUFixLogRatio = pufixlogratio
class JetSplitter_LowLogRatio (TrigJetSplitter):
__slots__ = []
......
......@@ -168,7 +168,7 @@ HLT::ErrorCode TrigBHremoval::hltExecute(std::vector<std::vector<HLT::TriggerEle
} else {
jetRatio = 999;
}
/*
if (jetEt < m_minJetEt) {
if (msgLvl() <= MSG::DEBUG)
msg() << MSG::DEBUG << "Jet "<< i << " below the " << m_minJetEt << " GeV threshold; Et " << jetEt << "; skipping this jet." << endreq;
......@@ -184,7 +184,7 @@ HLT::ErrorCode TrigBHremoval::hltExecute(std::vector<std::vector<HLT::TriggerEle
msg() << MSG::DEBUG << "Jet "<< i << " below the " << m_minLogRatio << " threshold for the log-ratio cut; logRatio = " << jetRatio << "; skipping this jet." << endreq;
continue;
}
*/
int countCaloCell=0;
// float min_dPhi=99;
//calculate number of cells in line of fire
......
......@@ -46,7 +46,8 @@ TrigJetSplitter::TrigJetSplitter(const std::string & name, ISvcLocator* pSvcLoca
declareProperty ("JetMinEt", m_minJetEt = 30.0); // in GeV (increase from 15 GeV to be same as vertex threshold)
declareProperty ("JetMaxEta", m_maxJetEta = 2.5+m_etaHalfWidth); // tracker acceptance + jet half-width
declareProperty ("JetLogRatio", m_logRatio = 1.2); // minimum/maximum log ratio for being selected
declareProperty("Reversed", m_reversedCut = false, "reversed cut for collimated photons");
declareProperty ("JetPUFixLogRatio", m_pufixLogRatio = 1.2); // minimum log ratio for pileup removal
declareProperty ("Reversed", m_reversedCut = false, "reversed cut for collimated photons");
}
......@@ -69,8 +70,9 @@ HLT::ErrorCode TrigJetSplitter::hltInitialize() {
msg() << MSG::DEBUG << " MinJetEt = " << m_minJetEt << endreq;
msg() << MSG::DEBUG << " MaxJetEta = " << m_maxJetEta << endreq;
msg() << MSG::DEBUG << " LogRatio = " << m_logRatio << endreq;
msg() << MSG::DEBUG << " PUFixLogRatio = " << m_pufixLogRatio << endreq;
msg() << MSG::DEBUG << " ReversedCut = " << m_reversedCut << endreq;
}
}
return HLT::OK;
}
......@@ -131,27 +133,27 @@ HLT::ErrorCode TrigJetSplitter::hltExecute(std::vector<std::vector<HLT::TriggerE
// Retrieve calo cells --> this is now done in TrigBHremoval.cxx
// -----------------------------------
/*
std::vector<const CaloCellContainer*> vectorOfCellContainers;
std::vector<const CaloCellContainer*> vectorOfCellContainers;
if(getFeatures(jetTE.front(), vectorOfCellContainers, "") != HLT::OK) {
if(getFeatures(jetTE.front(), vectorOfCellContainers, "") != HLT::OK) {
msg() << MSG::WARNING << "Failed to get TrigCells" << endreq;
return HLT::OK;
}
}
if(msgLvl() <= MSG::DEBUG) msg() << MSG::DEBUG << "Got vector with " << vectorOfCellContainers.size() << " CellContainers" << endreq;
if(msgLvl() <= MSG::DEBUG) msg() << MSG::DEBUG << "Got vector with " << vectorOfCellContainers.size() << " CellContainers" << endreq;
// if no containers were found, just leave the vector empty and leave
if ( vectorOfCellContainers.size() < 1) {
// if no containers were found, just leave the vector empty and leave
if ( vectorOfCellContainers.size() < 1) {
msg() << MSG::ERROR << "No cells to analyse, leaving!" << endreq;
return HLT::OK;
}
}
// get last container to be put in vector (should also be the only one)
const CaloCellContainer* theCellCont = vectorOfCellContainers.back();
// get last container to be put in vector (should also be the only one)
const CaloCellContainer* theCellCont = vectorOfCellContainers.back();
if(msgLvl() <= MSG::DEBUG) {
if(msgLvl() <= MSG::DEBUG) {
msg() << MSG::DEBUG << " Retrieved a Cell Container of Size= " << theCellCont->size() << endreq;
}
}
*/
if (msgLvl() <= MSG::DEBUG)
msg() << MSG::DEBUG << "Found " << jets->size() << " jets, creating corresponding RoIs" << endreq;
......@@ -188,12 +190,86 @@ HLT::ErrorCode TrigJetSplitter::hltExecute(std::vector<std::vector<HLT::TriggerE
continue;
}
if(!m_reversedCut) {
if ( jetRatio < m_logRatio) {
if ( jetRatio < m_pufixLogRatio) {
if (msgLvl() <= MSG::DEBUG)
msg() << MSG::DEBUG << "Jet "<< i << " below the " << m_logRatio << " threshold for the log-ratio cut; logRatio = " << jetRatio << "; skipping this jet." << endreq;
msg() << MSG::DEBUG << "Jet "<< i << " below the " << m_pufixLogRatio << " threshold for the log-ratio cut; logRatio = " << jetRatio << "; skipping this jet." << endreq;
continue;
}
} else {
double pufixLR = -1; // Will recalculate the jet's logRatio if it is not already > 1.2
if ( jetRatio < m_logRatio ) {
if (msgLvl() <= MSG::DEBUG)
msg() << MSG::DEBUG << "Jet "<< i << " is below the " << m_logRatio << " threshold for the log-ratio cut and above the " << m_pufixLogRatio << " threshold; logRatio = " << jetRatio << "; performing pileup removal." << endreq;
size_t nClusters = jet->numConstituents();
double clusterPU_sumEEM = 0; double clusterPU_sumE = 0;
for (size_t clust = 0; clust < nClusters; clust++) {
const xAOD::CaloCluster * aCluster = dynamic_cast<const xAOD::CaloCluster*> (jet->rawConstituent(clust));
double clusEEM = 0;
clusEEM+=(aCluster)->eSample(CaloSampling::EMB1);
clusEEM+=(aCluster)->eSample(CaloSampling::EMB2);
clusEEM+=(aCluster)->eSample(CaloSampling::EMB3);
clusEEM+=(aCluster)->eSample(CaloSampling::EME1);
clusEEM+=(aCluster)->eSample(CaloSampling::EME2);
clusEEM+=(aCluster)->eSample(CaloSampling::EME3);
clusEEM+=(aCluster)->eSample(CaloSampling::FCAL1);
double lambda = aCluster->getMomentValue(xAOD::CaloCluster::CENTER_LAMBDA);
if (lambda > 500) continue;
double d_eta = aCluster->rawEta() - jetEta;
double d_phi = HLT::wrapPhi(aCluster->rawPhi() - jetPhi);
double d_R2 = d_eta*d_eta + d_phi*d_phi;
if (d_R2 < 0.15*0.15) continue;
clusterPU_sumEEM+=clusEEM/1000.;
clusterPU_sumE+=aCluster->rawE()/1000.;
}
double jetEEM_EMscale = 0; double jetE_EMscale = 0; //Working on EM scale because calE() doesn't always return correct EEM and cluster moment EMF not accessable during testing
std::vector<double> samplingEnergy = jet->getAttribute<std::vector<double> >("EnergyPerSampling");
for(size_t s=0; s<samplingEnergy.size(); s++) {
double samplingE = 0.001*(samplingEnergy.at(s));
if ( s < 8 || (s > 20 && s < 28) ) jetEEM_EMscale+=samplingE; // EM layers 0-7 and 21-27
jetE_EMscale+=samplingE;
}
double pufixEMF = (jetEEM_EMscale - clusterPU_sumEEM)/(jetE_EMscale - clusterPU_sumE);
if (CxxUtils::fpcompare::greater(pufixEMF,zero)){
if(CxxUtils::fpcompare::greater_equal(pufixEMF,one)) pufixLR = -999.;
else pufixLR = log10(double(1./pufixEMF - 1.));
} else {
pufixLR = 999;
}
if ( pufixLR < m_logRatio) {
if (msgLvl() <= MSG::DEBUG)
msg() << MSG::DEBUG << "Jet "<< i << " is still below the " << m_logRatio << " threshold for the log-ratio cut; recalculated logRatio = " << pufixLR << "; skipping this jet." << endreq;
continue;
}
if (msgLvl() <= MSG::DEBUG){
msg() << MSG::DEBUG << "Jet "<< i << " old logRatio: " << jetRatio << " new logRatio: " << pufixLR << endreq;
msg() << MSG::DEBUG << "This is the jet with ET: " << jetEt << " eta: " << jetEta << " phi: " << jetPhi << endreq;
}
}
jetRatio = pufixLR;
}
else {
if ( jetRatio > m_logRatio) {
if (msgLvl() <= MSG::DEBUG)
msg() << MSG::DEBUG << "Jet "<< i << " above the " << m_logRatio << " threshold for the log-ratio cut; logRatio = " << jetRatio << "; skipping this jet." << endreq;
......@@ -208,34 +284,34 @@ HLT::ErrorCode TrigJetSplitter::hltExecute(std::vector<std::vector<HLT::TriggerE
//-tile
//-provenance
//-energy
if((*celliter)->caloDDE()->is_tile() && (*celliter)->energy() > 240.0){
if((*celliter)->caloDDE()->is_tile() && (*celliter)->energy() > 240.0){
//-dPhi
//-dR
float d_phi = HLT::wrapPhi((*celliter)->phi() - jetPhi);
float d_eta = (*celliter)->eta() - jetEta;
if(fabs(d_phi) < 0.2 && sqrt(d_eta*d_eta + d_phi*d_phi) > 0.3){
//-early
float t = (*celliter)->time();
if(t < -2.0){
//-halo time
float x = (*celliter)->x();
float y = (*celliter)->y();
float z = (*celliter)->z();
float c = 299.792458;//mm per ns
float r = sqrt(x*x + y*y);
if((fabs(t - (z-sqrt(z*z + r*r))/c) < 5.0) || (fabs(t - (-z-sqrt(z*z + r*r))/c) < 5.0)){
msg() << MSG::DEBUG << " cell E = " << (*celliter)->energy() << " cell phi = " << (*celliter)->phi() << " cell eta = " << (*celliter)->eta() << " cell r = " << r << endreq;
countCaloCell++;
}
}
}
}
//-dPhi
//-dR
float d_phi = HLT::wrapPhi((*celliter)->phi() - jetPhi);
float d_eta = (*celliter)->eta() - jetEta;
if(fabs(d_phi) < 0.2 && sqrt(d_eta*d_eta + d_phi*d_phi) > 0.3){
//-early
float t = (*celliter)->time();
if(t < -2.0){
//-halo time
float x = (*celliter)->x();
float y = (*celliter)->y();
float z = (*celliter)->z();
float c = 299.792458;//mm per ns
float r = sqrt(x*x + y*y);
if((fabs(t - (z-sqrt(z*z + r*r))/c) < 5.0) || (fabs(t - (-z-sqrt(z*z + r*r))/c) < 5.0)){
msg() << MSG::DEBUG << " cell E = " << (*celliter)->energy() << " cell phi = " << (*celliter)->phi() << " cell eta = " << (*celliter)->eta() << " cell r = " << r << endreq;
countCaloCell++;
}
}
}
}
}
if (msgLvl() <= MSG::DEBUG)
msg() << MSG::DEBUG << "Jet "<< i << "; Et " << jetEt << "; eta "<< jetEta << "; phi " << jetPhi <<"; logRatio " << jetRatio << "; LoF Cells " << countCaloCell << endreq;
msg() << MSG::DEBUG << "Jet "<< i << "; Et " << jetEt << "; eta "<< jetEta << "; phi " << jetPhi <<"; logRatio " << jetRatio << "; LoF Cells " << countCaloCell << endreq;
*/
// Create an output TE seeded by an empty vector
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment