| 1 |
#include "UserCode/yangyong/interface/EcalRecHitsRegional.h"
|
| 2 |
#include "DataFormats/EcalRecHit/interface/EcalRecHitCollections.h"
|
| 3 |
#include "DataFormats/HLTReco/interface/TriggerFilterObjectWithRefs.h"
|
| 4 |
#include "FWCore/MessageLogger/interface/MessageLogger.h"
|
| 5 |
#include "CondFormats/EcalObjects/interface/EcalChannelStatus.h"
|
| 6 |
#include "CondFormats/DataRecord/interface/EcalChannelStatusRcd.h"
|
| 7 |
#include "FWCore/Framework/interface/ESHandle.h"
|
| 8 |
|
| 9 |
EcalRecHitsRegional::EcalRecHitsRegional(const edm::ParameterSet& iConfig)
|
| 10 |
{
|
| 11 |
barrelHits_ = iConfig.getParameter< edm::InputTag > ("barrelHitCollection");
|
| 12 |
endcapHits_ = iConfig.getParameter< edm::InputTag > ("endcapHitCollection");
|
| 13 |
|
| 14 |
selectedBarrelHits_ =
|
| 15 |
iConfig.getParameter< std::string > ("selectedBarrelHitCollection");
|
| 16 |
selectedEndcapHits_ =
|
| 17 |
iConfig.getParameter< std::string > ("selectedEndcapHitCollection");
|
| 18 |
|
| 19 |
|
| 20 |
RegionalMatch_ = iConfig.getParameter<bool>("RegionalMatch");
|
| 21 |
ptMinEMObj_ = iConfig.getParameter<double>("ptMinEMObj");
|
| 22 |
EMregionEtaMargin_ = iConfig.getParameter<double>("EMregionEtaMargin"); ///0.25, //
|
| 23 |
EMregionPhiMargin_ = iConfig.getParameter<double>("EMregionPhiMargin"); ///0.4
|
| 24 |
l1IsolatedTag_ = iConfig.getParameter< edm::InputTag > ("l1IsolatedTag");
|
| 25 |
l1NonIsolatedTag_ = iConfig.getParameter< edm::InputTag > ("l1NonIsolatedTag");
|
| 26 |
debug_ = iConfig.getParameter<int> ("debugLevel");
|
| 27 |
|
| 28 |
|
| 29 |
TheMapping = new EcalElectronicsMapping();
|
| 30 |
first_ = true;
|
| 31 |
|
| 32 |
//register your products
|
| 33 |
produces< EBRecHitCollection >(selectedBarrelHits_);
|
| 34 |
produces< EERecHitCollection >(selectedEndcapHits_);
|
| 35 |
|
| 36 |
///produces<trigger::TriggerFilterObjectWithRefs>();
|
| 37 |
|
| 38 |
|
| 39 |
}
|
| 40 |
|
| 41 |
|
| 42 |
EcalRecHitsRegional::~EcalRecHitsRegional()
|
| 43 |
{
|
| 44 |
|
| 45 |
|
| 46 |
}
|
| 47 |
|
| 48 |
|
| 49 |
// ------------ method called to produce the data ------------
|
| 50 |
bool
|
| 51 |
EcalRecHitsRegional::filter(edm::Event& iEvent, const edm::EventSetup& iSetup)
|
| 52 |
{
|
| 53 |
|
| 54 |
|
| 55 |
using namespace edm;
|
| 56 |
using namespace std;
|
| 57 |
|
| 58 |
|
| 59 |
|
| 60 |
if (first_) {
|
| 61 |
edm::ESHandle< EcalElectronicsMapping > ecalmapping;
|
| 62 |
iSetup.get< EcalMappingRcd >().get(ecalmapping);
|
| 63 |
const EcalElectronicsMapping* TheMapping_ = ecalmapping.product();
|
| 64 |
*TheMapping = *TheMapping_;
|
| 65 |
first_ = false;
|
| 66 |
|
| 67 |
}
|
| 68 |
|
| 69 |
|
| 70 |
vector<int>::iterator it;
|
| 71 |
|
| 72 |
|
| 73 |
FEDListUsed.clear();
|
| 74 |
if( RegionalMatch_ ){
|
| 75 |
|
| 76 |
// Get the L1 CaloGeometry
|
| 77 |
edm::ESHandle<L1CaloGeometry> l1CaloGeom ;
|
| 78 |
iSetup.get<L1CaloGeometryRecord>().get(l1CaloGeom) ;
|
| 79 |
edm::Handle< l1extra::L1EmParticleCollection > l1EGIso;
|
| 80 |
edm::Handle< l1extra::L1EmParticleCollection > l1EGNonIso ;
|
| 81 |
///Get L1 EG ojbects
|
| 82 |
try{
|
| 83 |
iEvent.getByLabel(l1IsolatedTag_, l1EGIso ) ;
|
| 84 |
for( l1extra::L1EmParticleCollection::const_iterator emItr = l1EGIso->begin();
|
| 85 |
emItr != l1EGIso->end() ;++emItr ){
|
| 86 |
|
| 87 |
if( emItr -> pt() < ptMinEMObj_ ) continue;
|
| 88 |
|
| 89 |
|
| 90 |
if( debug_>=2) cout<<"L1EGISO: "<< emItr -> pt() <<" "<< emItr ->eta()<<" "<<emItr -> phi()<<endl;
|
| 91 |
|
| 92 |
int etaIndex = emItr->gctEmCand()->etaIndex() ;
|
| 93 |
int phiIndex = emItr->gctEmCand()->phiIndex() ;
|
| 94 |
double etaLow = l1CaloGeom->etaBinLowEdge( etaIndex ) ;
|
| 95 |
double etaHigh = l1CaloGeom->etaBinHighEdge( etaIndex ) ;
|
| 96 |
double phiLow = l1CaloGeom->emJetPhiBinLowEdge( phiIndex ) ;
|
| 97 |
double phiHigh = l1CaloGeom->emJetPhiBinHighEdge( phiIndex ) ;
|
| 98 |
|
| 99 |
std::vector<int> feds = ListOfFEDS(etaLow, etaHigh, phiLow, phiHigh, EMregionEtaMargin_, EMregionPhiMargin_);
|
| 100 |
for (int n=0; n < (int)feds.size(); n++) {
|
| 101 |
int fed = feds[n];
|
| 102 |
it = find(FEDListUsed.begin(),FEDListUsed.end(),fed);
|
| 103 |
if( it == FEDListUsed.end()){
|
| 104 |
FEDListUsed.push_back(fed);
|
| 105 |
}
|
| 106 |
}
|
| 107 |
}
|
| 108 |
|
| 109 |
|
| 110 |
}catch(std::exception& ex ){
|
| 111 |
cout<<"l1EGIso not working.."<<endl;
|
| 112 |
}
|
| 113 |
|
| 114 |
try{
|
| 115 |
iEvent.getByLabel(l1NonIsolatedTag_, l1EGNonIso ) ;
|
| 116 |
|
| 117 |
for( l1extra::L1EmParticleCollection::const_iterator emItr = l1EGNonIso->begin();
|
| 118 |
emItr != l1EGNonIso->end() ;++emItr ){
|
| 119 |
|
| 120 |
if( emItr -> pt()< ptMinEMObj_ ) continue;
|
| 121 |
|
| 122 |
if( debug_>=2) cout<<"L1EGnonIso: "<< emItr -> pt() <<" "<< emItr ->eta()<<" "<<emItr -> phi()<<endl;
|
| 123 |
|
| 124 |
int etaIndex = emItr->gctEmCand()->etaIndex() ;
|
| 125 |
int phiIndex = emItr->gctEmCand()->phiIndex() ;
|
| 126 |
double etaLow = l1CaloGeom->etaBinLowEdge( etaIndex ) ;
|
| 127 |
double etaHigh = l1CaloGeom->etaBinHighEdge( etaIndex ) ;
|
| 128 |
double phiLow = l1CaloGeom->emJetPhiBinLowEdge( phiIndex ) ;
|
| 129 |
double phiHigh = l1CaloGeom->emJetPhiBinHighEdge( phiIndex ) ;
|
| 130 |
|
| 131 |
std::vector<int> feds = ListOfFEDS(etaLow, etaHigh, phiLow, phiHigh, EMregionEtaMargin_, EMregionPhiMargin_);
|
| 132 |
for (int n=0; n < (int)feds.size(); n++) {
|
| 133 |
int fed = feds[n];
|
| 134 |
it = find(FEDListUsed.begin(),FEDListUsed.end(),fed);
|
| 135 |
if( it == FEDListUsed.end()){
|
| 136 |
FEDListUsed.push_back(fed);
|
| 137 |
}
|
| 138 |
}
|
| 139 |
}
|
| 140 |
|
| 141 |
}catch(std::exception& ex ){
|
| 142 |
cout<<"l1EGnonIso not working.."<<endl;
|
| 143 |
}
|
| 144 |
|
| 145 |
|
| 146 |
|
| 147 |
} //// end of getting FED List
|
| 148 |
///separate into barrel and endcap to speed up when checking
|
| 149 |
FEDListUsedBarrel.clear();
|
| 150 |
FEDListUsedEndcap.clear();
|
| 151 |
for( int j=0; j< int(FEDListUsed.size());j++){
|
| 152 |
int fed = FEDListUsed[j];
|
| 153 |
if( fed >= 10 && fed <= 45){
|
| 154 |
FEDListUsedBarrel.push_back(fed);
|
| 155 |
}else FEDListUsedEndcap.push_back(fed);
|
| 156 |
}
|
| 157 |
|
| 158 |
|
| 159 |
|
| 160 |
// edm::ESHandle<EcalChannelStatus> csHandle;
|
| 161 |
// if (! useRecoFlag_) iSetup.get<EcalChannelStatusRcd>().get(csHandle);
|
| 162 |
// const EcalChannelStatus& channelStatus = *csHandle;
|
| 163 |
|
| 164 |
|
| 165 |
|
| 166 |
Handle<EBRecHitCollection> barrelRecHitsHandle;
|
| 167 |
Handle<EERecHitCollection> endcapRecHitsHandle;
|
| 168 |
|
| 169 |
|
| 170 |
iEvent.getByLabel(barrelHits_,barrelRecHitsHandle);
|
| 171 |
iEvent.getByLabel(endcapHits_,endcapRecHitsHandle);
|
| 172 |
|
| 173 |
//Create empty output collections
|
| 174 |
std::auto_ptr< EBRecHitCollection > selectedEBRecHitCollection( new EBRecHitCollection );
|
| 175 |
std::auto_ptr< EERecHitCollection > selectedEERecHitCollection( new EERecHitCollection );
|
| 176 |
|
| 177 |
// The Filter object. We don't really need to put anything into it, but we
|
| 178 |
// write an empty one for consistency
|
| 179 |
std::auto_ptr<trigger::TriggerFilterObjectWithRefs>
|
| 180 |
filterproduct (new trigger::TriggerFilterObjectWithRefs(path(),module()));
|
| 181 |
|
| 182 |
|
| 183 |
//Select interesting EcalRecHits (barrel)
|
| 184 |
EBRecHitCollection::const_iterator itb;
|
| 185 |
for (itb=barrelRecHitsHandle->begin(); itb!=barrelRecHitsHandle->end(); itb++) {
|
| 186 |
|
| 187 |
|
| 188 |
EBDetId det = itb->id();
|
| 189 |
if (RegionalMatch_){
|
| 190 |
int fed = TheMapping->DCCid(det);
|
| 191 |
it = find(FEDListUsedBarrel.begin(),FEDListUsedBarrel.end(),fed);
|
| 192 |
if(it == FEDListUsedBarrel.end()) continue;
|
| 193 |
}
|
| 194 |
|
| 195 |
selectedEBRecHitCollection->push_back(*itb);
|
| 196 |
|
| 197 |
}
|
| 198 |
|
| 199 |
//Select interesting EcalRecHits (endcaps)
|
| 200 |
EERecHitCollection::const_iterator ite;
|
| 201 |
for (ite=endcapRecHitsHandle->begin(); ite!=endcapRecHitsHandle->end(); ite++) {
|
| 202 |
|
| 203 |
EEDetId det = ite->id();
|
| 204 |
|
| 205 |
if (RegionalMatch_){
|
| 206 |
EcalElectronicsId elid = TheMapping->getElectronicsId(det);
|
| 207 |
int fed = elid.dccId();
|
| 208 |
it = find(FEDListUsedEndcap.begin(),FEDListUsedEndcap.end(),fed);
|
| 209 |
if(it == FEDListUsedEndcap.end()) continue;
|
| 210 |
}
|
| 211 |
|
| 212 |
selectedEERecHitCollection->push_back(*ite);
|
| 213 |
|
| 214 |
}
|
| 215 |
|
| 216 |
if( debug_>=1) cout<<" selectedEB/EErechits: "<< selectedEBRecHitCollection->size()<<" "<<selectedEERecHitCollection->size()<<endl;
|
| 217 |
|
| 218 |
|
| 219 |
//Put selected information in the event
|
| 220 |
iEvent.put( selectedEBRecHitCollection, selectedBarrelHits_);
|
| 221 |
iEvent.put( selectedEERecHitCollection, selectedEndcapHits_);
|
| 222 |
|
| 223 |
return true;
|
| 224 |
|
| 225 |
}
|
| 226 |
|
| 227 |
|
| 228 |
|
| 229 |
//////FED list this is obsolete
|
| 230 |
std::vector<int> EcalRecHitsRegional::ListOfFEDS(double etaLow, double etaHigh, double phiLow,
|
| 231 |
double phiHigh, double etamargin, double phimargin)
|
| 232 |
{
|
| 233 |
|
| 234 |
std::vector<int> FEDs;
|
| 235 |
|
| 236 |
if (phimargin > Geom::pi()) phimargin = Geom::pi() ;
|
| 237 |
|
| 238 |
|
| 239 |
if (debug_>=2) std::cout << " etaLow etaHigh phiLow phiHigh " << etaLow << " " <<
|
| 240 |
etaHigh << " " << phiLow << " " << phiHigh << std::endl;
|
| 241 |
|
| 242 |
etaLow -= etamargin;
|
| 243 |
etaHigh += etamargin;
|
| 244 |
double phiMinus = phiLow - phimargin;
|
| 245 |
double phiPlus = phiHigh + phimargin;
|
| 246 |
|
| 247 |
bool all = false;
|
| 248 |
double dd = fabs(phiPlus-phiMinus);
|
| 249 |
if (debug_>=2) std::cout << " dd = " << dd << std::endl;
|
| 250 |
if (dd > 2.*Geom::pi() ) all = true;
|
| 251 |
|
| 252 |
while (phiPlus > Geom::pi()) { phiPlus -= 2.*Geom::pi() ; }
|
| 253 |
while (phiMinus < 0) { phiMinus += 2.*Geom::pi() ; }
|
| 254 |
if ( phiMinus > Geom::pi()) phiMinus -= 2.*Geom::pi() ;
|
| 255 |
|
| 256 |
double dphi = phiPlus - phiMinus;
|
| 257 |
if (dphi < 0) dphi += 2.*Geom::pi() ;
|
| 258 |
if (debug_>=2) std::cout << "dphi = " << dphi << std::endl;
|
| 259 |
if (dphi > Geom::pi()) {
|
| 260 |
int fed_low1 = TheMapping -> GetFED(etaLow,phiMinus*180./Geom::pi());
|
| 261 |
int fed_low2 = TheMapping -> GetFED(etaLow,phiPlus*180./Geom::pi());
|
| 262 |
if (debug_>=2) std::cout << "fed_low1 fed_low2 " << fed_low1 << " " << fed_low2 << std::endl;
|
| 263 |
if (fed_low1 == fed_low2) all = true;
|
| 264 |
int fed_hi1 = TheMapping -> GetFED(etaHigh,phiMinus*180./Geom::pi());
|
| 265 |
int fed_hi2 = TheMapping -> GetFED(etaHigh,phiPlus*180./Geom::pi());
|
| 266 |
if (debug_>=2) std::cout << "fed_hi1 fed_hi2 " << fed_hi1 << " " << fed_hi2 << std::endl;
|
| 267 |
if (fed_hi1 == fed_hi2) all = true;
|
| 268 |
}
|
| 269 |
|
| 270 |
if (all) {
|
| 271 |
if (debug_>=2) std::cout << " unpack everything in phi ! " << std::endl;
|
| 272 |
phiMinus = -20 * Geom::pi() / 180.; // -20 deg
|
| 273 |
phiPlus = -40 * Geom::pi() / 180.; // -20 deg
|
| 274 |
}
|
| 275 |
|
| 276 |
if (debug_>=2) std::cout << " with margins : " << etaLow << " " << etaHigh << " " <<
|
| 277 |
phiMinus << " " << phiPlus << std::endl;
|
| 278 |
|
| 279 |
|
| 280 |
const EcalEtaPhiRegion ecalregion(etaLow,etaHigh,phiMinus,phiPlus);
|
| 281 |
|
| 282 |
FEDs = TheMapping -> GetListofFEDs(ecalregion);
|
| 283 |
|
| 284 |
/*
|
| 285 |
if (debug_) {
|
| 286 |
int nn = (int)FEDs.size();
|
| 287 |
for (int ii=0; ii < nn; ii++) {
|
| 288 |
std::cout << "unpack fed " << FEDs[ii] << std::endl;
|
| 289 |
}
|
| 290 |
}
|
| 291 |
*/
|
| 292 |
|
| 293 |
return FEDs;
|
| 294 |
|
| 295 |
}
|
| 296 |
|
| 297 |
//define this as a plug-in
|
| 298 |
DEFINE_FWK_MODULE(EcalRecHitsRegional);
|
