00001 #include "RecoTauTag/RecoTau/interface/TCTauAlgorithm.h"
00002
00003 #include "DataFormats/TrackReco/interface/Track.h"
00004 #include "DataFormats/VertexReco/interface/Vertex.h"
00005 #include "RecoTauTag/TauTagTools/interface/TauTagTools.h"
00006 #include "Geometry/Records/interface/CaloGeometryRecord.h"
00007
00008 #include "TrackingTools/TransientTrack/interface/TransientTrackBuilder.h"
00009
00010 #include "Math/VectorUtil.h"
00011 using namespace ROOT::Math;
00012
00013 using namespace reco;
00014
00015 TCTauAlgorithm::TCTauAlgorithm(){
00016 init();
00017 }
00018
00019 TCTauAlgorithm::TCTauAlgorithm(const edm::ParameterSet& iConfig){
00020 init();
00021 inputConfig(iConfig);
00022 }
00023
00024 TCTauAlgorithm::~TCTauAlgorithm(){}
00025
00026 void TCTauAlgorithm::init(){
00027
00028 event = 0;
00029 setup = 0;
00030
00031 trackAssociator = new TrackDetectorAssociator();
00032 trackAssociator->useDefaultPropagator();
00033
00034 all = 0;
00035 passed = 0;
00036 prongs = -1;
00037 algoComponentUsed = 0;
00038 }
00039
00040 void TCTauAlgorithm::inputConfig(const edm::ParameterSet& iConfig){
00041
00042 etCaloOverTrackMin = iConfig.getParameter<double>("EtCaloOverTrackMin");
00043 etCaloOverTrackMax = iConfig.getParameter<double>("EtCaloOverTrackMax");
00044 etHcalOverTrackMin = iConfig.getParameter<double>("EtHcalOverTrackMin");
00045 etHcalOverTrackMax = iConfig.getParameter<double>("EtHcalOverTrackMax");
00046
00047 signalCone = iConfig.getParameter<double>("SignalConeSize");
00048 ecalCone = iConfig.getParameter<double>("EcalConeSize");
00049
00050 EcalRecHitsEB_input= iConfig.getParameter<edm::InputTag>("EBRecHitCollection");
00051 EcalRecHitsEE_input= iConfig.getParameter<edm::InputTag>("EERecHitCollection");
00052 HBHERecHits_input = iConfig.getParameter<edm::InputTag>("HBHERecHitCollection");
00053 HORecHits_input = iConfig.getParameter<edm::InputTag>("HORecHitCollection");
00054 HFRecHits_input = iConfig.getParameter<edm::InputTag>("HFRecHitCollection");
00055
00056 edm::ParameterSet pset = iConfig.getParameter<edm::ParameterSet>("TrackAssociatorParameters");
00057 trackAssociatorParameters.loadParameters( pset );
00058
00059 dropCaloJets = iConfig.getUntrackedParameter<bool>("DropCaloJets",false);
00060 dropRejected = iConfig.getUntrackedParameter<bool>("DropRejectedJets",true);
00061 }
00062
00063
00064 double TCTauAlgorithm::efficiency(){
00065 return double(passed)/all;
00066 }
00067
00068 int TCTauAlgorithm::statistics(){
00069 return passed;
00070 }
00071
00072 int TCTauAlgorithm::allTauCandidates(){
00073 return all;
00074 }
00075
00076 int TCTauAlgorithm::algoComponent(){
00077 return algoComponentUsed;
00078 }
00079
00080 void TCTauAlgorithm::eventSetup(const edm::Event& iEvent,const edm::EventSetup& iSetup){
00081
00082 event = &iEvent;
00083 setup = &iSetup;
00084
00085 edm::ESHandle<TransientTrackBuilder> builder;
00086 iSetup.get<TransientTrackRecord>().get("TransientTrackBuilder",builder);
00087 transientTrackBuilder = builder.product();
00088
00089
00090 edm::ESHandle<CaloGeometry> geometry;
00091 iSetup.get<CaloGeometryRecord>().get(geometry);
00092
00093
00094 EB = geometry->getSubdetectorGeometry(DetId::Ecal,EcalBarrel);
00095 EE = geometry->getSubdetectorGeometry(DetId::Ecal,EcalEndcap);
00096 HB = geometry->getSubdetectorGeometry(DetId::Hcal,HcalBarrel);
00097 HE = geometry->getSubdetectorGeometry(DetId::Hcal,HcalEndcap);
00098 HO = geometry->getSubdetectorGeometry(DetId::Hcal,HcalOuter);
00099 HF = geometry->getSubdetectorGeometry(DetId::Hcal,HcalForward);
00100
00101
00102 iEvent.getByLabel( EcalRecHitsEB_input, EBRecHits );
00103 iEvent.getByLabel( EcalRecHitsEE_input, EERecHits );
00104
00105 iEvent.getByLabel( HBHERecHits_input, HBHERecHits );
00106 iEvent.getByLabel( HORecHits_input, HORecHits );
00107 iEvent.getByLabel( HFRecHits_input, HFRecHits );
00108 }
00109
00110
00111 math::XYZTLorentzVector TCTauAlgorithm::recalculateEnergy(const reco::CaloTau& jet){
00112
00113 const TrackRef& leadTk = jet.leadTrack();
00114
00115 const TrackRefVector associatedTracks = jet.caloTauTagInfoRef()->Tracks();
00116
00117 const CaloJet* cJet = jet.caloTauTagInfoRef()->calojetRef().get();
00118 CaloJet caloJet = *cJet;
00119 caloJet.setP4(jet.p4());
00120
00121 return recalculateEnergy(caloJet,leadTk,associatedTracks);
00122 }
00123
00124 math::XYZTLorentzVector TCTauAlgorithm::recalculateEnergy(const reco::CaloJet& caloJet,const TrackRef& leadTk,const TrackRefVector& associatedTracks){
00125
00126 all++;
00127
00128 math::XYZTLorentzVector p4(0,0,0,0);
00129 algoComponentUsed = TCAlgoUndetermined;
00130
00131
00132 p4 = caloJet.p4();
00133
00134 if(leadTk.isNull()) return p4;
00135
00136 XYZVector momentum(0,0,0);
00137 int prongCounter = 0;
00138 edm::RefVector<TrackCollection>::const_iterator iTrack;
00139 for(iTrack = associatedTracks.begin(); iTrack!= associatedTracks.end(); ++iTrack){
00140 double DR = ROOT::Math::VectorUtil::DeltaR(leadTk->momentum(),(*iTrack)->momentum());
00141 if(DR < signalCone) {
00142 momentum+=(*iTrack)->momentum();
00143 prongCounter++;
00144 }
00145 }
00146 if(momentum.Rho() == 0) return p4;
00147
00148 XYZVector ltrackEcalHitPoint = trackEcalHitPoint(*leadTk);
00149
00150 if(! (ltrackEcalHitPoint.Rho() > 0 && ltrackEcalHitPoint.Rho() < 9999) ) return p4;
00151
00152 std::pair<XYZVector,XYZVector> caloClusters = getClusterEnergy(caloJet,ltrackEcalHitPoint,signalCone);
00153 XYZVector EcalCluster = caloClusters.first;
00154 XYZVector HcalCluster = caloClusters.second;
00155
00156 double eCaloOverTrack = (EcalCluster.R()+HcalCluster.R()-momentum.R())/momentum.R();
00157
00158 std::pair<XYZVector,XYZVector> caloClustersPhoton = getClusterEnergy(caloJet,ltrackEcalHitPoint,ecalCone);
00159 XYZVector EcalClusterPhoton = caloClustersPhoton.first;
00160
00161 math::XYZTLorentzVector p4photons(0,0,0,EcalClusterPhoton.R() - EcalCluster.R());
00162
00163 if( eCaloOverTrack > etCaloOverTrackMin && eCaloOverTrack < etCaloOverTrackMax ) {
00164
00165 double eHcalOverTrack = (HcalCluster.R()-momentum.R())/momentum.R();
00166
00167 if ( eHcalOverTrack > etHcalOverTrackMin && eHcalOverTrack < etHcalOverTrackMax ) {
00168 p4.SetXYZT(EcalCluster.X() + momentum.X(),
00169 EcalCluster.Y() + momentum.Y(),
00170 EcalCluster.Z() + momentum.Z(),
00171 EcalCluster.R() + momentum.R());
00172 p4 += p4photons;
00173 algoComponentUsed = TCAlgoMomentumECAL;
00174 }else{
00175 p4.SetXYZT(momentum.X(),
00176 momentum.Y(),
00177 momentum.Z(),
00178 momentum.R());
00179 algoComponentUsed = TCAlgoMomentum;
00180 }
00181 }
00182 if( eCaloOverTrack > etCaloOverTrackMax ) {
00183 double eHcalOverTrack = (HcalCluster.R()-momentum.R())/momentum.R();
00184
00185 if ( eHcalOverTrack > etHcalOverTrackMin && eHcalOverTrack < etHcalOverTrackMax ) {
00186 p4.SetXYZT(EcalCluster.X() + momentum.X(),
00187 EcalCluster.Y() + momentum.Y(),
00188 EcalCluster.Z() + momentum.Z(),
00189 EcalCluster.R() + momentum.R());
00190 p4 += p4photons;
00191 algoComponentUsed = TCAlgoMomentumECAL;
00192 }
00193 if ( eHcalOverTrack < etHcalOverTrackMin ) {
00194 if(!dropCaloJets) p4.SetXYZT(caloJet.px(),caloJet.py(),caloJet.pz(),caloJet.energy());
00195 else p4.SetXYZT(0,0,0,0);
00196 algoComponentUsed = TCAlgoCaloJet;
00197 }
00198 if ( eHcalOverTrack > etHcalOverTrackMax ) {
00199 algoComponentUsed = TCAlgoHadronicJet;
00200 if(!dropRejected) p4.SetXYZT(caloJet.px(),caloJet.py(),caloJet.pz(),caloJet.energy());
00201 else p4.SetXYZT(0,0,0,0);
00202 }
00203 }
00204 if( eCaloOverTrack < etCaloOverTrackMin ) {
00205 algoComponentUsed = TCAlgoTrackProblem;
00206 if(!dropRejected) p4.SetXYZT(caloJet.px(),caloJet.py(),caloJet.pz(),caloJet.energy());
00207 }
00208
00209 if(p4.Et() > 0) passed++;
00210
00211 return p4;
00212 }
00213
00214
00215 XYZVector TCTauAlgorithm::trackEcalHitPoint(const TransientTrack& transientTrack,const CaloJet& caloJet){
00216
00217
00218 GlobalPoint ecalHitPosition(0,0,0);
00219
00220 double maxTowerEt = 0;
00221 std::vector<CaloTowerPtr> towers = caloJet.getCaloConstituents();
00222 for(std::vector<CaloTowerPtr>::const_iterator iTower = towers.begin();
00223 iTower!= towers.end(); ++iTower){
00224 if((*iTower)->et() > maxTowerEt){
00225 maxTowerEt = (*iTower)->et();
00226 ecalHitPosition = (*iTower)->emPosition();
00227 }
00228 }
00229
00230
00231 XYZVector ecalHitPoint(0,0,0);
00232
00233 GlobalPoint trackEcalHitPoint = transientTrack.stateOnSurface(ecalHitPosition).globalPosition();
00234
00235 ecalHitPoint.SetXYZ(trackEcalHitPoint.x(),
00236 trackEcalHitPoint.y(),
00237 trackEcalHitPoint.z());
00238
00239 return ecalHitPoint;
00240 }
00241
00242 XYZVector TCTauAlgorithm::trackEcalHitPoint(const Track& track){
00243
00244 const FreeTrajectoryState fts = trackAssociator->getFreeTrajectoryState(*setup,track);
00245 TrackDetMatchInfo info = trackAssociator->associate(*event, *setup, fts, trackAssociatorParameters);
00246 if( info.isGoodEcal != 0 ) {
00247 return XYZVector(info.trkGlobPosAtEcal.x(),info.trkGlobPosAtEcal.y(),info.trkGlobPosAtEcal.z());
00248 }
00249 return XYZVector(0,0,0);
00250 }
00251
00252 std::pair<XYZVector,XYZVector> TCTauAlgorithm::getClusterEnergy(const CaloJet& caloJet,XYZVector& trackEcalHitPoint,double cone){
00253
00254 XYZVector ecalCluster(0,0,0);
00255 XYZVector hcalCluster(0,0,0);
00256
00257 std::vector<CaloTowerPtr> towers = caloJet.getCaloConstituents();
00258
00259 for(std::vector<CaloTowerPtr>::const_iterator iTower = towers.begin();
00260 iTower!= towers.end(); ++iTower){
00261 std::vector<XYZVector> ECALCells;
00262 std::vector<XYZVector> HCALCells;
00263
00264 size_t numRecHits = (**iTower).constituentsSize();
00265
00266
00267 for(size_t j = 0; j < numRecHits; j++) {
00268 DetId recHitDetID = (**iTower).constituent(j);
00269
00270 if( recHitDetID.det() == DetId::Ecal ){
00271 if( recHitDetID.subdetId() == 1 ){
00272 EBDetId ecalID = recHitDetID;
00273 EBRecHitCollection::const_iterator theRecHit = EBRecHits->find(ecalID);
00274 if(theRecHit != EBRecHits->end()){
00275 DetId id = theRecHit->detid();
00276 const CaloCellGeometry* this_cell = EB->getGeometry(id);
00277 double energy = theRecHit->energy();
00278 ECALCells.push_back(getCellMomentum(this_cell,energy));
00279 }
00280 }
00281 if( recHitDetID.subdetId() == 2 ){
00282 EEDetId ecalID = recHitDetID;
00283 EERecHitCollection::const_iterator theRecHit = EERecHits->find(ecalID);
00284 if(theRecHit != EERecHits->end()){
00285 DetId id = theRecHit->detid();
00286 const CaloCellGeometry* this_cell = EE->getGeometry(id);
00287 double energy = theRecHit->energy();
00288 ECALCells.push_back(getCellMomentum(this_cell,energy));
00289 }
00290 }
00291 }
00292 if( recHitDetID.det() == DetId::Hcal ){
00293 HcalDetId hcalID = recHitDetID;
00294 if( recHitDetID.subdetId() == HcalBarrel ){
00295
00296
00297 HBHERecHitCollection::const_iterator theRecHit=HBHERecHits->find(hcalID);
00298 if(theRecHit != HBHERecHits->end()){
00299 DetId id = theRecHit->detid();
00300 const CaloCellGeometry* this_cell = HB->getGeometry(id);
00301 double energy = theRecHit->energy();
00302 HCALCells.push_back(getCellMomentum(this_cell,energy));
00303 }
00304
00305 }
00306 if( recHitDetID.subdetId() == HcalEndcap ){
00307
00308
00309 HBHERecHitCollection::const_iterator theRecHit=HBHERecHits->find(hcalID);
00310 if(theRecHit != HBHERecHits->end()){
00311 DetId id = theRecHit->detid();
00312 const CaloCellGeometry* this_cell = HE->getGeometry(id);
00313 double energy = theRecHit->energy();
00314 HCALCells.push_back(getCellMomentum(this_cell,energy));
00315 }
00316
00317 }
00318 if( recHitDetID.subdetId() == HcalOuter ){
00319 HORecHitCollection::const_iterator theRecHit=HORecHits->find(hcalID);
00320 if(theRecHit != HORecHits->end()){
00321 DetId id = theRecHit->detid();
00322 const CaloCellGeometry* this_cell = HO->getGeometry(id);
00323 double energy = theRecHit->energy();
00324 HCALCells.push_back(getCellMomentum(this_cell,energy));
00325 }
00326 }
00327 if( recHitDetID.subdetId() == HcalForward ){
00328 HFRecHitCollection::const_iterator theRecHit=HFRecHits->find(hcalID);
00329 if(theRecHit != HFRecHits->end()){
00330 DetId id = theRecHit->detid();
00331 const CaloCellGeometry* this_cell = HF->getGeometry(id);
00332 double energy = theRecHit->energy();
00333 HCALCells.push_back(getCellMomentum(this_cell,energy));
00334 }
00335 }
00336 }
00337 }
00338
00339 std::vector<XYZVector>::const_iterator i;
00340 for(i = ECALCells.begin(); i != ECALCells.end(); ++i) {
00341 double DR = ROOT::Math::VectorUtil::DeltaR(trackEcalHitPoint,*i);
00342 if( DR < cone ) ecalCluster += *i;
00343 }
00344 for(i = HCALCells.begin(); i != HCALCells.end(); ++i) {
00345 double DR = ROOT::Math::VectorUtil::DeltaR(trackEcalHitPoint,*i);
00346 if( DR < cone ) hcalCluster += *i;
00347 }
00348 }
00349 return std::pair<XYZVector,XYZVector> (ecalCluster,hcalCluster);
00350 }
00351
00352 XYZVector TCTauAlgorithm::getCellMomentum(const CaloCellGeometry* cell,double& energy){
00353 XYZVector momentum(0,0,0);
00354 if(cell){
00355 GlobalPoint hitPosition = cell->getPosition();
00356
00357 double phi = hitPosition.phi();
00358 double theta = hitPosition.theta();
00359 if(theta > 3.14159) theta = 2*3.14159 - theta;
00360 double px = energy * sin(theta)*cos(phi);
00361 double py = energy * sin(theta)*sin(phi);
00362 double pz = energy * cos(theta);
00363
00364 momentum = XYZVector(px,py,pz);
00365 }
00366 return momentum;
00367 }
1.7.3