| 1 |
#include "../interface/PhotonIsolator.h"
|
| 2 |
|
| 3 |
PhotonIsolator::PhotonIsolator(edm::ParameterSet * config, edm::ParameterSet * producersNames):config_(config), producersNames_(producersNames)
|
| 4 |
{
|
| 5 |
verbosity_ = (*config_).getUntrackedParameter<int>("verbosity", 0);
|
| 6 |
basicClustersIsolation_BarrelBC_type_ = (*config_).getParameter<int>("basicClustersIsolation_BarrelBC_type");
|
| 7 |
basicClustersIsolation_EndcapBC_type_ = (*config_).getParameter<int>("basicClustersIsolation_EndcapBC_type");
|
| 8 |
basicClustersIsolation_DRmax_ = (*config_).getParameter<double>("basicClustersIsolation_DRmax");
|
| 9 |
basicClustersIsolation_ClusterEt_threshold_ = (*config_).getParameter<double>("basicClustersIsolation_ClusterEt_threshold");
|
| 10 |
basicClustersDoubleConeIsolation_BarrelBC_type_ = (*config_).getParameter<int>("basicClustersDoubleConeIsolation_BarrelBC_type");
|
| 11 |
basicClustersDoubleConeIsolation_EndcapBC_type_ = (*config_).getParameter<int>("basicClustersDoubleConeIsolation_EndcapBC_type");
|
| 12 |
basicClustersDoubleConeIsolation_DRmin_ = (*config_).getParameter<double>("basicClustersDoubleConeIsolation_DRmin");
|
| 13 |
basicClustersDoubleConeIsolation_DRmax_ = (*config_).getParameter<double>("basicClustersDoubleConeIsolation_DRmax");
|
| 14 |
basicClustersDoubleConeIsolation_ClusterEt_threshold_ = (*config_).getParameter<double>("basicClustersDoubleConeIsolation_ClusterEt_threshold");
|
| 15 |
hcalRecHitIsolation_DRmax_ = (*config_).getParameter<double>("hcalRecHitIsolation_DRmax");
|
| 16 |
hcalRecHitIsolation_HitEt_threshold_ = (*config_).getParameter<double>("hcalRecHitIsolation_HitEt_threshold");
|
| 17 |
trackerIsolation_DRmax_ = (*config_).getParameter<double>("trackerIsolation_DRmax");
|
| 18 |
trackerIsolation_pt_threshold_ = (*config_).getParameter<double>("trackerIsolation_pt_threshold");
|
| 19 |
trackerIsolation_pixelLayers_threshold_ = (*config_).getParameter<int>("trackerIsolation_pixelLayers_threshold");
|
| 20 |
allowMissingCollection_ = (*producersNames_).getUntrackedParameter<bool>("allowMissingCollection", false);
|
| 21 |
}
|
| 22 |
|
| 23 |
|
| 24 |
PhotonIsolator::~PhotonIsolator()
|
| 25 |
{
|
| 26 |
}
|
| 27 |
|
| 28 |
|
| 29 |
float PhotonIsolator::deltaR(float phi1, float phi2, float eta1, float eta2)
|
| 30 |
{
|
| 31 |
float dphi=deltaPhi(phi1,phi2);
|
| 32 |
float deta=fabs(eta1-eta2);
|
| 33 |
float dr = sqrt(dphi*dphi+ deta*deta);
|
| 34 |
return dr;
|
| 35 |
}
|
| 36 |
|
| 37 |
|
| 38 |
float PhotonIsolator::deltaPhi(float phi1, float phi2)
|
| 39 |
{
|
| 40 |
float dphi;
|
| 41 |
if(phi1<0) phi1+=TWOPI;
|
| 42 |
if(phi2<0) phi2+=TWOPI;
|
| 43 |
dphi=fabs(phi1-phi2);
|
| 44 |
if(dphi>TWOPI) dphi-=TWOPI;
|
| 45 |
if(dphi>PI) dphi=TWOPI-dphi;
|
| 46 |
return dphi;
|
| 47 |
}
|
| 48 |
|
| 49 |
|
| 50 |
bool PhotonIsolator::loadHcalRecHits(const edm::Event& iEvent, const edm::EventSetup& iSetup)
|
| 51 |
{
|
| 52 |
// Get HB/HE rechits
|
| 53 |
try
|
| 54 |
{
|
| 55 |
edm::InputTag hbheRecHitProducer = producersNames_->getParameter<edm::InputTag>("hbheRecHitProducer");
|
| 56 |
edm::Handle<HBHERecHitCollection> hbheHandle;
|
| 57 |
iEvent.getByLabel(hbheRecHitProducer, hbheHandle);
|
| 58 |
hbheHits_ = hbheHandle.product();
|
| 59 |
}
|
| 60 |
catch (cms::Exception& exception)
|
| 61 |
{
|
| 62 |
if ( !allowMissingCollection_ )
|
| 63 |
{
|
| 64 |
cout << " ##### ERROR IN PhotonIsolator::LoadHcalRecHits => No HBHERecHitCollection #####"<<endl;
|
| 65 |
throw exception;
|
| 66 |
}
|
| 67 |
if(verbosity_>1) cout << " ===> No HBHERecHitCollection, skip HCAL rechits isolation" << endl;
|
| 68 |
return false;
|
| 69 |
}
|
| 70 |
|
| 71 |
// Get HO rechits
|
| 72 |
try
|
| 73 |
{
|
| 74 |
edm::InputTag hoRecHitProducer = producersNames_->getParameter<edm::InputTag>("hoRecHitProducer");
|
| 75 |
edm::Handle<HORecHitCollection> hoHandle;
|
| 76 |
iEvent.getByLabel(hoRecHitProducer, hoHandle);
|
| 77 |
hoHits_ = hoHandle.product();
|
| 78 |
}
|
| 79 |
catch (cms::Exception& exception)
|
| 80 |
{
|
| 81 |
if ( !allowMissingCollection_ )
|
| 82 |
{
|
| 83 |
cout << " ##### ERROR IN PhotonIsolator::LoadHcalRecHits => No HORecHitCollection #####"<<endl;
|
| 84 |
throw exception;
|
| 85 |
}
|
| 86 |
if(verbosity_>1) cout << " ===> No HORecHitCollection, skip HCAL rechits isolation" << endl;
|
| 87 |
return false;
|
| 88 |
}
|
| 89 |
|
| 90 |
// Get HF rechits
|
| 91 |
try
|
| 92 |
{
|
| 93 |
edm::InputTag hfRecHitProducer = producersNames_->getParameter<edm::InputTag>("hfRecHitProducer");
|
| 94 |
edm::Handle<HFRecHitCollection> hfHandle;
|
| 95 |
iEvent.getByLabel(hfRecHitProducer, hfHandle);
|
| 96 |
hfHits_ = hfHandle.product();
|
| 97 |
}
|
| 98 |
catch (cms::Exception& exception)
|
| 99 |
{
|
| 100 |
if ( !allowMissingCollection_ )
|
| 101 |
{
|
| 102 |
cout << " ##### ERROR IN PhotonIsolator::LoadHcalRecHits => No HFRecHitCollection #####"<<endl;
|
| 103 |
throw exception;
|
| 104 |
}
|
| 105 |
if(verbosity_>1) cout << " ===> No HFRecHitCollection, skip HCAL rechits isolation" << endl;
|
| 106 |
return false;
|
| 107 |
}
|
| 108 |
|
| 109 |
// Get Calo Geometry
|
| 110 |
try
|
| 111 |
{
|
| 112 |
edm::ESHandle<CaloGeometry> geomHandle;
|
| 113 |
(iSetup.get<CaloGeometryRecord>()).get(geomHandle);
|
| 114 |
geometry_ = geomHandle.product();
|
| 115 |
}
|
| 116 |
catch (cms::Exception& exception)
|
| 117 |
{
|
| 118 |
if ( !allowMissingCollection_ )
|
| 119 |
{
|
| 120 |
cout << " ##### ERROR IN PhotonIsolator::LoadHcalRecHits => No CaloGeometryRecord #####"<<endl;
|
| 121 |
throw exception;
|
| 122 |
}
|
| 123 |
if(verbosity_>1) cout << " ===> No CaloGeometryRecord, skip HCAL rechits isolation" << endl;
|
| 124 |
return false;
|
| 125 |
}
|
| 126 |
|
| 127 |
return true;
|
| 128 |
}
|
| 129 |
|
| 130 |
|
| 131 |
Double_t PhotonIsolator::basicClustersIsolation(TRootPhoton* photon, TClonesArray* superClusters, TClonesArray* basicClusters)
|
| 132 |
{
|
| 133 |
Double_t sumEt = 0.;
|
| 134 |
Int_t isc_barrel = photon->scIndexOfType(basicClustersIsolation_BarrelBC_type_);
|
| 135 |
Int_t isc_endcap = photon->scIndexOfType(basicClustersIsolation_EndcapBC_type_);
|
| 136 |
Int_t isc = max(isc_barrel,isc_endcap);
|
| 137 |
|
| 138 |
if(isc==-1) { cout << " ##### ERROR IN PhotonIsolator::BasicClustersIsolation => no supercluster associated to the photon #####" << endl; return -1.;}
|
| 139 |
|
| 140 |
TRootSuperCluster* photonSC = (TRootSuperCluster*) superClusters->At(isc);
|
| 141 |
//cout << "Associated SC[" << isc << "]" << endl;;
|
| 142 |
|
| 143 |
TRootCluster* localBC;
|
| 144 |
for (int ibc=0; ibc<basicClusters->GetEntriesFast(); ibc++)
|
| 145 |
{
|
| 146 |
localBC = (TRootCluster*) basicClusters->At(ibc);
|
| 147 |
//cout << "localBC[" << ibc << "] type=" << localBC->type() << endl;
|
| 148 |
if( localBC->type()==basicClustersIsolation_BarrelBC_type_ || localBC->type()==basicClustersIsolation_EndcapBC_type_ )
|
| 149 |
{
|
| 150 |
if( localBC->Pt()>basicClustersIsolation_ClusterEt_threshold_ )
|
| 151 |
{
|
| 152 |
// No vertex correction for BC, so used uncorrected SC position for isolation
|
| 153 |
//Double_t dR = deltaR( photon->Phi(), localBC->Phi(), photon->Eta(), localBC->Eta() );
|
| 154 |
Double_t dR = deltaR( photonSC->Phi(), localBC->Phi(), photonSC->Eta(), localBC->Eta() );
|
| 155 |
//cout << "localBC DeltaR=" << dR << endl;
|
| 156 |
if( dR<basicClustersIsolation_DRmax_ )
|
| 157 |
{
|
| 158 |
bool inSuperCluster = false;
|
| 159 |
for (unsigned int isub=0; isub<photonSC->subBasicClusterIndexVector().size(); isub++)
|
| 160 |
{
|
| 161 |
// Veto Basic Clusters belonging to photon supercluster
|
| 162 |
if ( ibc==photonSC->subBasicClusterIndexVector().at(isub) )
|
| 163 |
{
|
| 164 |
inSuperCluster = true;
|
| 165 |
//break; //FIXME - Add break
|
| 166 |
}
|
| 167 |
}
|
| 168 |
//if (inSuperCluster) cout << "localBC in photonSC" << endl;
|
| 169 |
//if (!inSuperCluster) cout << "localBC not in photonSC" << endl;
|
| 170 |
if (!inSuperCluster) sumEt+=localBC->Pt();
|
| 171 |
//cout << "sumEt=" << sumEt << endl;
|
| 172 |
}
|
| 173 |
}
|
| 174 |
}
|
| 175 |
}
|
| 176 |
return sumEt;
|
| 177 |
}
|
| 178 |
|
| 179 |
|
| 180 |
Double_t PhotonIsolator::basicClustersDoubleConeIsolation(TRootPhoton* photon, TClonesArray* superClusters, TClonesArray* basicClusters)
|
| 181 |
{
|
| 182 |
Int_t isc_barrel = photon->scIndexOfType(basicClustersDoubleConeIsolation_BarrelBC_type_);
|
| 183 |
Int_t isc_endcap = photon->scIndexOfType(basicClustersDoubleConeIsolation_EndcapBC_type_);
|
| 184 |
Int_t isc = max(isc_barrel,isc_endcap);
|
| 185 |
|
| 186 |
if(isc==-1) { cout << " ##### ERROR IN PhotonIsolator::BasicClustersDoubleConeIsolation => no supercluster associated to the photon #####" << endl; return -1.;}
|
| 187 |
|
| 188 |
TRootSuperCluster* photonSC = (TRootSuperCluster*) superClusters->At(isc);
|
| 189 |
|
| 190 |
Double_t sumEt = 0.;
|
| 191 |
TRootCluster* localBC;
|
| 192 |
for (int ibc=0; ibc<basicClusters->GetEntriesFast(); ibc++)
|
| 193 |
{
|
| 194 |
localBC = (TRootCluster*) basicClusters->At(ibc);
|
| 195 |
if( localBC->type()==basicClustersDoubleConeIsolation_BarrelBC_type_ || localBC->type()==basicClustersDoubleConeIsolation_EndcapBC_type_ )
|
| 196 |
{
|
| 197 |
if( localBC->Pt()>basicClustersDoubleConeIsolation_ClusterEt_threshold_ )
|
| 198 |
{
|
| 199 |
// No vertex correction for BC, so used uncorrected SC position for isolation
|
| 200 |
//Double_t dR = deltaR( photon->Phi(), localBC->Phi(), photon->Eta(), localBC->Eta() );
|
| 201 |
Double_t dR = deltaR( photonSC->Phi(), localBC->Phi(), photonSC->Eta(), localBC->Eta() );
|
| 202 |
if( dR<basicClustersDoubleConeIsolation_DRmax_ && dR>basicClustersDoubleConeIsolation_DRmin_ )
|
| 203 |
{
|
| 204 |
sumEt+=localBC->Pt();
|
| 205 |
}
|
| 206 |
}
|
| 207 |
}
|
| 208 |
}
|
| 209 |
return sumEt;
|
| 210 |
}
|
| 211 |
|
| 212 |
|
| 213 |
Double_t PhotonIsolator::hcalRecHitIsolation(TRootPhoton* photon)
|
| 214 |
{
|
| 215 |
Double_t sumEt=0.;
|
| 216 |
|
| 217 |
// No vertex correction for recHits, so used uncorrected photon position
|
| 218 |
//double phi=photon->Phi();
|
| 219 |
//double eta=photon->Eta();
|
| 220 |
double phi=photon->caloPosition().Phi();
|
| 221 |
double eta=photon->caloPosition().Eta();
|
| 222 |
|
| 223 |
for(HBHERecHitCollection::const_iterator hbheItr = hbheHits_->begin(); hbheItr != hbheHits_->end(); hbheItr++)
|
| 224 |
{
|
| 225 |
double HcalHit_energy=hbheItr->energy();
|
| 226 |
double HcalHit_eta=geometry_->getPosition(hbheItr->id()).eta();
|
| 227 |
double HcalHit_phi=geometry_->getPosition(hbheItr->id()).phi();
|
| 228 |
float HcalHit_pth=HcalHit_energy*sin(2*atan(exp(-HcalHit_eta)));
|
| 229 |
|
| 230 |
if(HcalHit_pth>hcalRecHitIsolation_HitEt_threshold_)
|
| 231 |
{
|
| 232 |
float newDelta=deltaR(HcalHit_phi,phi,HcalHit_eta,eta);
|
| 233 |
if(newDelta<hcalRecHitIsolation_DRmax_) sumEt+=HcalHit_pth;
|
| 234 |
}
|
| 235 |
}
|
| 236 |
|
| 237 |
|
| 238 |
for(HFRecHitCollection::const_iterator hfItr = hfHits_->begin(); hfItr != hfHits_->end(); hfItr++)
|
| 239 |
{
|
| 240 |
double HcalHit_energy=hfItr->energy();
|
| 241 |
double HcalHit_eta=geometry_->getPosition(hfItr->id()).eta();
|
| 242 |
double HcalHit_phi=geometry_->getPosition(hfItr->id()).phi();
|
| 243 |
float HcalHit_pth=HcalHit_energy*sin(2*atan(exp(-HcalHit_eta)));
|
| 244 |
|
| 245 |
if(HcalHit_pth>hcalRecHitIsolation_HitEt_threshold_)
|
| 246 |
{
|
| 247 |
float newDelta=deltaR(HcalHit_phi,phi,HcalHit_eta,eta);
|
| 248 |
if(newDelta<hcalRecHitIsolation_DRmax_) sumEt+=HcalHit_pth;
|
| 249 |
}
|
| 250 |
}
|
| 251 |
|
| 252 |
|
| 253 |
for(HORecHitCollection::const_iterator hoItr = hoHits_->begin(); hoItr != hoHits_->end(); hoItr++)
|
| 254 |
{
|
| 255 |
double HcalHit_energy=hoItr->energy();
|
| 256 |
double HcalHit_eta=geometry_->getPosition(hoItr->id()).eta();
|
| 257 |
double HcalHit_phi=geometry_->getPosition(hoItr->id()).phi();
|
| 258 |
float HcalHit_pth=HcalHit_energy*sin(2*atan(exp(-HcalHit_eta)));
|
| 259 |
if(HcalHit_pth>hcalRecHitIsolation_HitEt_threshold_)
|
| 260 |
{
|
| 261 |
float newDelta=deltaR(HcalHit_phi,phi,HcalHit_eta,eta);
|
| 262 |
if(newDelta<hcalRecHitIsolation_DRmax_) sumEt+=HcalHit_pth;
|
| 263 |
}
|
| 264 |
}
|
| 265 |
|
| 266 |
return sumEt;
|
| 267 |
}
|
| 268 |
|
| 269 |
|
| 270 |
Double_t PhotonIsolator::trackerIsolation(TRootPhoton* photon, TClonesArray* tracks, Int_t &nTracks)
|
| 271 |
{
|
| 272 |
// FIXME - If photon not vertex corrected then correct it
|
| 273 |
Double_t sumPt = 0.;
|
| 274 |
nTracks = 0;
|
| 275 |
TRootTrack* localTrack;
|
| 276 |
for (int itk=0; itk<tracks->GetEntriesFast(); itk++)
|
| 277 |
{
|
| 278 |
localTrack = (TRootTrack*) tracks->At(itk);
|
| 279 |
if( localTrack->Pt() > trackerIsolation_pt_threshold_ && localTrack->pixelLayersWithMeasurement() >= trackerIsolation_pixelLayers_threshold_ )
|
| 280 |
{
|
| 281 |
Double_t dR = deltaR( photon->Phi(), localTrack->Phi(), photon->Eta(), localTrack->Eta() );
|
| 282 |
if( dR<trackerIsolation_DRmax_ )
|
| 283 |
{
|
| 284 |
sumPt+=localTrack->Pt();
|
| 285 |
nTracks++;
|
| 286 |
}
|
| 287 |
}
|
| 288 |
}
|
| 289 |
return sumPt;
|
| 290 |
}
|
