00001 #include "JetMETCorrections/TauJet/interface/HardTauAlgorithm.h"
00002
00003 #include "DataFormats/TrackReco/interface/Track.h"
00004 #include "DataFormats/VertexReco/interface/Vertex.h"
00005 #include "RecoTauTag/TauTagTools/interface/TauTagTools.h"
00006
00007 #include "TrackingTools/TransientTrack/interface/TransientTrackBuilder.h"
00008
00009 #include "Math/VectorUtil.h"
00010 using namespace ROOT::Math;
00011
00012 HardTauAlgorithm::HardTauAlgorithm(){
00013 init();
00014 }
00015
00016 HardTauAlgorithm::HardTauAlgorithm(const edm::ParameterSet& iConfig){
00017 init();
00018 inputConfig(iConfig);
00019 }
00020
00021 HardTauAlgorithm::~HardTauAlgorithm(){}
00022
00023 void HardTauAlgorithm::init(){
00024
00025 etCaloOverTrackMin = -0.9;
00026 etCaloOverTrackMax = 0.0;
00027 etHcalOverTrackMin = -0.5;
00028 etHcalOverTrackMax = 0.5;
00029
00030 signalCone = 0.2;
00031 ecalCone = 0.5;
00032 matchingCone = 0.1;
00033 tkptmin = 1.0;
00034
00035 tkmaxipt = 0.03;
00036 tkmaxChi2 = 100;
00037 tkminPixelHitsn = 2;
00038 tkminTrackerHitsn = 8;
00039
00040 trackInput = InputTag("generalTracks");
00041 vertexInput = InputTag("offlinePrimaryVertices");
00042
00043 EcalRecHitsEB_input = InputTag("ecalRecHit:EcalRecHitsEB");
00044 EcalRecHitsEE_input = InputTag("ecalRecHit:EcalRecHitsEE");
00045 HBHERecHits_input = InputTag("hbhereco");
00046 HORecHits_input = InputTag("horeco");
00047 HFRecHits_input = InputTag("hfreco");
00048
00049 all = 0;
00050 passed = 0;
00051 prongs = -1;
00052 }
00053
00054 void HardTauAlgorithm::inputConfig(const edm::ParameterSet& iConfig){
00055
00056 etCaloOverTrackMin = iConfig.getUntrackedParameter<double>("EtCaloOverTrackMin",etCaloOverTrackMin);
00057 etCaloOverTrackMax = iConfig.getUntrackedParameter<double>("EtCaloOverTrackMax",etCaloOverTrackMax);
00058 etHcalOverTrackMin = iConfig.getUntrackedParameter<double>("EtHcalOverTrackMin",etHcalOverTrackMin);
00059 etHcalOverTrackMax = iConfig.getUntrackedParameter<double>("EtHcalOverTrackMax",etHcalOverTrackMax);
00060
00061 signalCone = iConfig.getUntrackedParameter<double>("SignalConeSize",signalCone);
00062 ecalCone = iConfig.getUntrackedParameter<double>("EcalConeSize",ecalCone);
00063 matchingCone = iConfig.getUntrackedParameter<double>("MatchingConeSize",matchingCone);
00064 tkptmin = iConfig.getUntrackedParameter<double>("Track_minPt",tkptmin);
00065
00066 tkmaxipt = iConfig.getUntrackedParameter<double>("tkmaxipt",tkmaxipt);
00067 tkmaxChi2 = iConfig.getUntrackedParameter<double>("tkmaxChi2",tkmaxChi2);
00068 tkminPixelHitsn = iConfig.getUntrackedParameter<int>("tkminPixelHitsn",tkminPixelHitsn);
00069 tkminTrackerHitsn = iConfig.getUntrackedParameter<int>("tkminTrackerHitsn",tkminTrackerHitsn);
00070
00071 trackInput = iConfig.getUntrackedParameter<InputTag>("TrackCollection",trackInput);
00072 vertexInput = iConfig.getUntrackedParameter<InputTag>("PVProducer",vertexInput);
00073
00074 EcalRecHitsEB_input= iConfig.getUntrackedParameter<InputTag>("EBRecHitCollection",EcalRecHitsEB_input);
00075 EcalRecHitsEE_input= iConfig.getUntrackedParameter<InputTag>("EERecHitCollection",EcalRecHitsEE_input);
00076 HBHERecHits_input = iConfig.getUntrackedParameter<InputTag>("HBHERecHitCollection",HBHERecHits_input);
00077 HORecHits_input = iConfig.getUntrackedParameter<InputTag>("HORecHitCollection",HORecHits_input);
00078 HFRecHits_input = iConfig.getUntrackedParameter<InputTag>("HFRecHitCollection",HFRecHits_input);
00079
00080 }
00081
00082
00083 double HardTauAlgorithm::efficiency(){
00084 return double(passed)/all;
00085 }
00086
00087
00088 void HardTauAlgorithm::eventSetup(const edm::Event& iEvent,const edm::EventSetup& iSetup){
00089 edm::ESHandle<TransientTrackBuilder> builder;
00090 iSetup.get<TransientTrackRecord>().get("TransientTrackBuilder",builder);
00091 transientTrackBuilder = builder.product();
00092
00093
00094 ESHandle<CaloGeometry> geometry;
00095 iSetup.get<IdealGeometryRecord>().get(geometry);
00096
00097 EB = geometry->getSubdetectorGeometry(DetId::Ecal,EcalBarrel);
00098 EE = geometry->getSubdetectorGeometry(DetId::Ecal,EcalEndcap);
00099 HB = geometry->getSubdetectorGeometry(DetId::Hcal,HcalBarrel);
00100 HE = geometry->getSubdetectorGeometry(DetId::Hcal,HcalEndcap);
00101 HO = geometry->getSubdetectorGeometry(DetId::Hcal,HcalOuter);
00102 HF = geometry->getSubdetectorGeometry(DetId::Hcal,HcalForward);
00103
00104
00105 iEvent.getByLabel( EcalRecHitsEB_input, EBRecHits );
00106 iEvent.getByLabel( EcalRecHitsEE_input, EERecHits );
00107
00108 iEvent.getByLabel( HBHERecHits_input, HBHERecHits );
00109 iEvent.getByLabel( HORecHits_input, HORecHits );
00110 iEvent.getByLabel( HFRecHits_input, HFRecHits );
00111
00112
00113 iEvent.getByLabel(trackInput,tracks);
00114 iEvent.getByLabel(vertexInput,thePVs);
00115 }
00116
00117
00118 TLorentzVector HardTauAlgorithm::recalculateEnergy(const reco::CaloTau& jet){
00119
00120 const TrackRef& leadTk = jet.leadTrack();
00121
00122 const TrackRefVector associatedTracks = jet.caloTauTagInfoRef()->Tracks();
00123
00124 const CaloJet* cJet = jet.caloTauTagInfoRef()->calojetRef().get();
00125 CaloJet caloJet = *cJet;
00126 caloJet.setP4(jet.p4());
00127
00128 return recalculateEnergy(caloJet,leadTk,associatedTracks);
00129 }
00130
00131 TLorentzVector HardTauAlgorithm::recalculateEnergy(const reco::CaloJet& caloJet){
00132
00133 TrackRef leadTk;
00134 TrackRefVector associatedTracks;
00135
00136 double ptmax = 0;
00137
00138 if(tracks.isValid()) {
00139
00140
00141 TrackCollection::const_iterator iTrack;
00142 for(unsigned int i = 0; i < tracks->size(); ++i){
00143 TrackRef trackRef(tracks,i);
00144 double DR = ROOT::Math::VectorUtil::DeltaR(caloJet.momentum(),trackRef->momentum());
00145 if(DR < 0.5) associatedTracks.push_back(trackRef);
00146 }
00147 }
00148
00149 Vertex thePV = *(thePVs->begin());
00150 TrackRefVector theFilteredTracks = TauTagTools::filteredTracks(associatedTracks,
00151 tkptmin,
00152 tkminPixelHitsn,
00153 tkminTrackerHitsn,
00154 tkmaxipt,
00155 tkmaxChi2,
00156 thePV);
00157
00158 for(TrackRefVector::const_iterator i = theFilteredTracks.begin();
00159 i!= theFilteredTracks.end(); ++i){
00160 double DR = ROOT::Math::VectorUtil::DeltaR(caloJet.momentum(),(*i)->momentum());
00161 if(DR < matchingCone && (*i)->pt() > ptmax){
00162 leadTk = *i;
00163 ptmax = (*i)->pt();
00164 }
00165 }
00166
00167 if(ptmax > 0) return recalculateEnergy(caloJet,leadTk,theFilteredTracks);
00168
00169 return TLorentzVector(0,0,0,0);
00170 }
00171
00172 TLorentzVector HardTauAlgorithm::recalculateEnergy(const reco::Jet& tau){
00173
00174 cout << "HardTauAlgorithm::recalculateEnergy(const reco::Jet&) "
00175 << "is not working. " << endl;
00176 cout << "Please use CaloJet or CaloTau instead. Exiting..." << endl;
00177 exit(0);
00178
00179 const CaloJet& cJet = dynamic_cast<const CaloJet&>(tau);
00180 CaloJet caloJet = cJet;
00181 caloJet.setP4(tau.p4());
00182
00183 return recalculateEnergy(caloJet);
00184 }
00185
00186 TLorentzVector HardTauAlgorithm::recalculateEnergy(const reco::IsolatedTauTagInfo& tau){
00187
00188 const TrackRef& leadTk = tau.leadingSignalTrack(matchingCone,tkptmin);
00189
00190 const TrackRefVector associatedTracks = tau.allTracks();
00191
00192 const CaloJet& cJet = dynamic_cast<const CaloJet&>(*(tau.jet()));
00193 CaloJet caloJet = cJet;
00194 caloJet.setP4(tau.jet().get()->p4());
00195
00196 return recalculateEnergy(caloJet,leadTk,associatedTracks);
00197 }
00198
00199 TLorentzVector HardTauAlgorithm::recalculateEnergy(const reco::CaloJet& caloJet,const TrackRef& leadTk,const TrackRefVector& associatedTracks){
00200
00201 all++;
00202
00203 TLorentzVector p4(0,0,0,0);
00204
00205 if(leadTk->pt() == 0) return p4;
00206
00207 XYZVector momentum(0,0,0);
00208 int prongCounter = 0;
00209 RefVector<TrackCollection>::const_iterator iTrack;
00210 for(iTrack = associatedTracks.begin(); iTrack!= associatedTracks.end(); ++iTrack){
00211 double DR = ROOT::Math::VectorUtil::DeltaR(leadTk->momentum(),(*iTrack)->momentum());
00212 if(DR < signalCone) {
00213 momentum+=(*iTrack)->momentum();
00214 prongCounter++;
00215 }
00216 }
00217 if(momentum.Rho() == 0) return p4;
00218
00219 const TransientTrack transientTrack = transientTrackBuilder->build(*leadTk);
00220 XYZVector ltrackEcalHitPoint = trackEcalHitPoint(transientTrack,caloJet);
00221 if(! (ltrackEcalHitPoint.Rho() > 0 && ltrackEcalHitPoint.Rho() < 9999) ) return p4;
00222
00223 pair<XYZVector,XYZVector> caloClusters = getClusterEnergy(caloJet,ltrackEcalHitPoint,signalCone);
00224 XYZVector EcalCluster = caloClusters.first;
00225 XYZVector HcalCluster = caloClusters.second;
00226
00227 double etCaloOverTrack = (EcalCluster.Rho()+HcalCluster.Rho()-momentum.Rho())/momentum.Rho();
00228
00229 pair<XYZVector,XYZVector> caloClustersPhoton = getClusterEnergy(caloJet,ltrackEcalHitPoint,ecalCone);
00230 XYZVector EcalClusterPhoton = caloClustersPhoton.first;
00231 TLorentzVector p4photons(EcalClusterPhoton.X() - EcalCluster.X(),
00232 EcalClusterPhoton.Y() - EcalCluster.Y(),
00233 EcalClusterPhoton.Z() - EcalCluster.Z(),
00234 EcalClusterPhoton.R() - EcalCluster.R());
00235
00236 if( etCaloOverTrack > etCaloOverTrackMin && etCaloOverTrack < etCaloOverTrackMax ) {
00237 p4.SetXYZT(momentum.X(),
00238 momentum.Y(),
00239 momentum.Z(),
00240 momentum.R());
00241 p4 += p4photons;
00242 }
00243 if( etCaloOverTrack > etCaloOverTrackMax ) {
00244 double etHcalOverTrack = (HcalCluster.Rho()-momentum.Rho())/momentum.Rho();
00245
00246 if ( etHcalOverTrack > etHcalOverTrackMin && etHcalOverTrack < etHcalOverTrackMax ) {
00247 p4.SetXYZT(EcalCluster.X() + momentum.X(),
00248 EcalCluster.Y() + momentum.Y(),
00249 EcalCluster.Z() + momentum.Z(),
00250 EcalCluster.R() + momentum.R());
00251 p4 += p4photons;
00252 }
00253 if ( etHcalOverTrack < etHcalOverTrackMin ) {
00254 p4.SetXYZT(caloJet.px(),caloJet.py(),caloJet.pz(),caloJet.energy());
00255 }
00256 }
00257
00258 if(p4.Et() > 0) passed++;
00259
00260 return p4;
00261 }
00262
00263
00264 XYZVector HardTauAlgorithm::trackEcalHitPoint(const TransientTrack& transientTrack,const CaloJet& caloJet){
00265
00266
00267 GlobalPoint ecalHitPosition(0,0,0);
00268
00269 double maxTowerEt = 0;
00270 vector<CaloTowerPtr> towers = caloJet.getCaloConstituents();
00271 for(vector<CaloTowerPtr>::const_iterator iTower = towers.begin();
00272 iTower!= towers.end(); ++iTower){
00273 if((*iTower)->et() > maxTowerEt){
00274 maxTowerEt = (*iTower)->et();
00275 ecalHitPosition = GlobalPoint((*iTower)->momentum().x(),
00276 (*iTower)->momentum().y(),
00277 (*iTower)->momentum().z());
00278 }
00279 }
00280
00281
00282 XYZVector ecalHitPoint(0,0,0);
00283
00284 try{
00285 TrajectoryStateClosestToPoint TSCP = transientTrack.trajectoryStateClosestToPoint(ecalHitPosition);
00286 GlobalPoint trackEcalHitPoint = TSCP.position();
00287
00288 ecalHitPoint.SetXYZ(trackEcalHitPoint.x(),
00289 trackEcalHitPoint.y(),
00290 trackEcalHitPoint.z());
00291 }catch(...) {;}
00292
00293 return ecalHitPoint;
00294 }
00295
00296 pair<XYZVector,XYZVector> HardTauAlgorithm::getClusterEnergy(const CaloJet& caloJet,XYZVector& trackEcalHitPoint,double cone){
00297
00298 XYZVector ecalCluster(0,0,0);
00299 XYZVector hcalCluster(0,0,0);
00300
00301 vector<CaloTowerPtr> towers = caloJet.getCaloConstituents();
00302
00303 for(vector<CaloTowerPtr>::const_iterator iTower = towers.begin();
00304 iTower!= towers.end(); ++iTower){
00305 vector<XYZVector> ECALCells;
00306 vector<XYZVector> HCALCells;
00307
00308 size_t numRecHits = (**iTower).constituentsSize();
00309
00310
00311 for(size_t j = 0; j < numRecHits; j++) {
00312 DetId recHitDetID = (**iTower).constituent(j);
00313
00314 if( recHitDetID.det() == DetId::Ecal ){
00315 if( recHitDetID.subdetId() == 1 ){
00316 EBDetId ecalID = recHitDetID;
00317 EBRecHitCollection::const_iterator theRecHit = EBRecHits->find(ecalID);
00318 if(theRecHit != EBRecHits->end()){
00319 DetId id = theRecHit->detid();
00320 const CaloCellGeometry* this_cell = EB->getGeometry(id);
00321 double energy = theRecHit->energy();
00322 ECALCells.push_back(getCellMomentum(this_cell,energy));
00323 }
00324 }
00325 if( recHitDetID.subdetId() == 2 ){
00326 EEDetId ecalID = recHitDetID;
00327 EERecHitCollection::const_iterator theRecHit = EERecHits->find(ecalID);
00328 if(theRecHit != EERecHits->end()){
00329 DetId id = theRecHit->detid();
00330 const CaloCellGeometry* this_cell = EE->getGeometry(id);
00331 double energy = theRecHit->energy();
00332 ECALCells.push_back(getCellMomentum(this_cell,energy));
00333 }
00334 }
00335 }
00336 if( recHitDetID.det() == DetId::Hcal ){
00337 HcalDetId hcalID = recHitDetID;
00338 if( recHitDetID.subdetId() == HcalBarrel ){
00339
00340
00341 HBHERecHitCollection::const_iterator theRecHit=HBHERecHits->find(hcalID);
00342 if(theRecHit != HBHERecHits->end()){
00343 DetId id = theRecHit->detid();
00344 const CaloCellGeometry* this_cell = HB->getGeometry(id);
00345 double energy = theRecHit->energy();
00346 HCALCells.push_back(getCellMomentum(this_cell,energy));
00347 }
00348
00349 }
00350 if( recHitDetID.subdetId() == HcalEndcap ){
00351
00352
00353 HBHERecHitCollection::const_iterator theRecHit=HBHERecHits->find(hcalID);
00354 if(theRecHit != HBHERecHits->end()){
00355 DetId id = theRecHit->detid();
00356 const CaloCellGeometry* this_cell = HE->getGeometry(id);
00357 double energy = theRecHit->energy();
00358 HCALCells.push_back(getCellMomentum(this_cell,energy));
00359 }
00360
00361 }
00362 if( recHitDetID.subdetId() == HcalOuter ){
00363 HORecHitCollection::const_iterator theRecHit=HORecHits->find(hcalID);
00364 if(theRecHit != HORecHits->end()){
00365 DetId id = theRecHit->detid();
00366 const CaloCellGeometry* this_cell = HO->getGeometry(id);
00367 double energy = theRecHit->energy();
00368 HCALCells.push_back(getCellMomentum(this_cell,energy));
00369 }
00370 }
00371 if( recHitDetID.subdetId() == HcalForward ){
00372 HFRecHitCollection::const_iterator theRecHit=HFRecHits->find(hcalID);
00373 if(theRecHit != HFRecHits->end()){
00374 DetId id = theRecHit->detid();
00375 const CaloCellGeometry* this_cell = HF->getGeometry(id);
00376 double energy = theRecHit->energy();
00377 HCALCells.push_back(getCellMomentum(this_cell,energy));
00378 }
00379 }
00380 }
00381 }
00382
00383 vector<XYZVector>::const_iterator i;
00384 for(i = ECALCells.begin(); i != ECALCells.end(); ++i) {
00385 double DR = ROOT::Math::VectorUtil::DeltaR(trackEcalHitPoint,*i);
00386 if( DR < cone ) ecalCluster += *i;
00387 }
00388 for(i = HCALCells.begin(); i != HCALCells.end(); ++i) {
00389 double DR = ROOT::Math::VectorUtil::DeltaR(trackEcalHitPoint,*i);
00390 if( DR < cone ) hcalCluster += *i;
00391 }
00392 }
00393 return pair<XYZVector,XYZVector> (ecalCluster,hcalCluster);
00394 }
00395
00396 XYZVector HardTauAlgorithm::getCellMomentum(const CaloCellGeometry* cell,double& energy){
00397 XYZVector momentum(0,0,0);
00398 if(cell){
00399 GlobalPoint hitPosition = cell->getPosition();
00400
00401 double phi = hitPosition.phi();
00402 double theta = hitPosition.theta();
00403 if(theta > 3.14159) theta = 2*3.14159 - theta;
00404 double px = energy * sin(theta)*cos(phi);
00405 double py = energy * sin(theta)*sin(phi);
00406 double pz = energy * cos(theta);
00407
00408 momentum = XYZVector(px,py,pz);
00409 }
00410 return momentum;
00411 }
1.5.4