| 1 |
lethuill |
1.4 |
#include "../interface/PhotonIsolator.h"
|
| 2 |
mlethuil |
1.1 |
|
| 3 |
lethuill |
1.8 |
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 |
mlethuil |
1.1 |
|
| 50 |
lethuill |
1.8 |
bool PhotonIsolator::loadHcalRecHits(const edm::Event& iEvent, const edm::EventSetup& iSetup)
|
| 51 |
mlethuil |
1.1 |
{
|
| 52 |
lethuill |
1.8 |
// 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 |
mlethuil |
1.2 |
}
|
| 129 |
lethuill |
1.6 |
|
| 130 |
lethuill |
1.8 |
|
| 131 |
|
|
Double_t PhotonIsolator::basicClustersIsolation(TRootPhoton* photon, TClonesArray* superClusters, TClonesArray* basicClusters)
|
| 132 |
|
|
{
|
| 133 |
|
|
Double_t sumEt = 0.;
|
| 134 |
lethuill |
1.9 |
Int_t isc_barrel = photon->scIndexOfType(basicClustersIsolation_BarrelBC_type_);
|
| 135 |
|
|
Int_t isc_endcap = photon->scIndexOfType(basicClustersIsolation_EndcapBC_type_);
|
| 136 |
lethuill |
1.8 |
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 |
lethuill |
1.9 |
for (unsigned int isub=0; isub<photonSC->subBasicClusterIndexVector().size(); isub++)
|
| 160 |
lethuill |
1.8 |
{
|
| 161 |
|
|
// Veto Basic Clusters belonging to photon supercluster
|
| 162 |
lethuill |
1.9 |
if ( ibc==photonSC->subBasicClusterIndexVector().at(isub) )
|
| 163 |
lethuill |
1.8 |
{
|
| 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 |
lethuill |
1.9 |
Int_t isc_barrel = photon->scIndexOfType(basicClustersDoubleConeIsolation_BarrelBC_type_);
|
| 183 |
|
|
Int_t isc_endcap = photon->scIndexOfType(basicClustersDoubleConeIsolation_EndcapBC_type_);
|
| 184 |
lethuill |
1.8 |
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 |
|
|
}
|
