CMS 3D CMS Logo

CMSSW_4_4_3_patch1/src/DQMOffline/Trigger/src/EgHLTOffHelper.cc

Go to the documentation of this file.
00001 #include "DQMOffline/Trigger/interface/EgHLTOffHelper.h"
00002 
00003 
00004 #include "DataFormats/EcalDetId/interface/EcalSubdetector.h"
00005 #include "DataFormats/DetId/interface/DetId.h"
00006 #include "DataFormats/JetReco/interface/CaloJetCollection.h"
00007 #include "FWCore/Common/interface/TriggerNames.h"
00008 
00009 #include "RecoEcal/EgammaCoreTools/interface/EcalClusterTools.h"
00010 
00011 #include "DataFormats/TrackReco/interface/Track.h"
00012 #include "DataFormats/GsfTrackReco/interface/GsfTrack.h"
00013 
00014 #include "RecoEgamma/EgammaHLTAlgos/interface/EgammaHLTTrackIsolation.h"
00015 #include "RecoEgamma/EgammaIsolationAlgos/interface/EgammaTowerIsolation.h"
00016 #include "RecoEgamma/EgammaIsolationAlgos/interface/EgammaRecHitIsolation.h"
00017 #include "RecoLocalCalo/EcalRecAlgos/interface/EcalSeverityLevelAlgo.h"
00018 #include "RecoLocalCalo/EcalRecAlgos/interface/EcalSeverityLevelAlgoRcd.h"
00019 
00020 #include "DQMOffline/Trigger/interface/EgHLTTrigCodes.h"
00021 #include "DQMOffline/Trigger/interface/EgHLTTrigTools.h"
00022 #include "DQMOffline/Trigger/interface/EgHLTErrCodes.h"
00023 
00024 #include <iostream>
00025 
00026 using namespace egHLT;
00027 
00028 
00029 
00030 OffHelper::~OffHelper()
00031 {
00032   if(hltEleTrkIsolAlgo_) delete hltEleTrkIsolAlgo_;
00033   if(hltPhoTrkIsolAlgo_) delete hltPhoTrkIsolAlgo_;
00034 }
00035 
00036 void OffHelper::setup(const edm::ParameterSet& conf)
00037 {
00038 
00039   ecalRecHitsEBTag_ = conf.getParameter<edm::InputTag>("BarrelRecHitCollection");
00040   ecalRecHitsEETag_ = conf.getParameter<edm::InputTag>("EndcapRecHitCollection");
00041   caloJetsTag_ = conf.getParameter<edm::InputTag>("CaloJetCollection");
00042   isolTrkTag_ = conf.getParameter<edm::InputTag>("IsolTrackCollection");
00043   hbheHitsTag_ = conf.getParameter<edm::InputTag>("HBHERecHitCollection");
00044   hfHitsTag_ = conf.getParameter<edm::InputTag>("HFRecHitCollection");
00045   electronsTag_ = conf.getParameter<edm::InputTag>("ElectronCollection");
00046   photonsTag_ = conf.getParameter<edm::InputTag>("PhotonCollection");
00047   triggerSummaryLabel_ = conf.getParameter<edm::InputTag>("triggerSummaryLabel");
00048   hltTag_ = conf.getParameter<std::string>("hltTag");
00049   beamSpotTag_ = conf.getParameter<edm::InputTag>("BeamSpotProducer");
00050   caloTowersTag_ = conf.getParameter<edm::InputTag>("CaloTowers");
00051   trigResultsTag_ = conf.getParameter<edm::InputTag>("TrigResults");
00052   vertexTag_ = conf.getParameter<edm::InputTag>("VertexCollection");
00053 
00054   eleCuts_.setup(conf.getParameter<edm::ParameterSet>("eleCuts"));
00055   eleLooseCuts_.setup(conf.getParameter<edm::ParameterSet>("eleLooseCuts"));
00056   phoCuts_.setup(conf.getParameter<edm::ParameterSet>("phoCuts"));
00057   phoLooseCuts_.setup(conf.getParameter<edm::ParameterSet>("phoLooseCuts"));
00058  
00059   //now we have the isolations completely configurable via python
00060   hltEMIsolOuterCone_ = conf.getParameter<double>("hltEMIsolOuterCone");
00061   hltEMIsolInnerConeEB_ = conf.getParameter<double>("hltEMIsolInnerConeEB");
00062   hltEMIsolEtaSliceEB_ = conf.getParameter<double>("hltEMIsolEtaSliceEB");
00063   hltEMIsolEtMinEB_ = conf.getParameter<double>("hltEMIsolEtMinEB");
00064   hltEMIsolEMinEB_ = conf.getParameter<double>("hltEMIsolEMinEB");
00065   hltEMIsolInnerConeEE_ = conf.getParameter<double>("hltEMIsolInnerConeEE");
00066   hltEMIsolEtaSliceEE_ = conf.getParameter<double>("hltEMIsolEtaSliceEE");
00067   hltEMIsolEtMinEE_ = conf.getParameter<double>("hltEMIsolEtMinEE");
00068   hltEMIsolEMinEE_ = conf.getParameter<double>("hltEMIsolEMinEE");
00069 
00070   hltPhoTrkIsolPtMin_ = conf.getParameter<double>("hltPhoTrkIsolPtMin");
00071   hltPhoTrkIsolOuterCone_ = conf.getParameter<double>("hltPhoTrkIsolOuterCone");
00072   hltPhoTrkIsolInnerCone_ = conf.getParameter<double>("hltPhoTrkIsolInnerCone"); 
00073   hltPhoTrkIsolZSpan_ = conf.getParameter<double>("hltPhoTrkIsolZSpan");
00074   hltPhoTrkIsolRSpan_ = conf.getParameter<double>("hltPhoTrkIsolZSpan");
00075   hltPhoTrkIsolCountTrks_ = conf.getParameter<bool>("hltPhoTrkIsolCountTrks");
00076 
00077   hltEleTrkIsolPtMin_ = conf.getParameter<double>("hltEleTrkIsolPtMin");
00078   hltEleTrkIsolOuterCone_ = conf.getParameter<double>("hltEleTrkIsolOuterCone");
00079   hltEleTrkIsolInnerCone_ = conf.getParameter<double>("hltEleTrkIsolInnerCone"); 
00080   hltEleTrkIsolZSpan_ = conf.getParameter<double>("hltEleTrkIsolZSpan");
00081   hltEleTrkIsolRSpan_ = conf.getParameter<double>("hltEleTrkIsolZSpan");
00082 
00083   hltHadIsolOuterCone_ = conf.getParameter<double>("hltHadIsolOuterCone");
00084   hltHadIsolInnerCone_ = conf.getParameter<double>("hltHadIsolInnerCone");
00085   hltHadIsolEtMin_ = conf.getParameter<double>("hltHadIsolEtMin");
00086   hltHadIsolDepth_ = conf.getParameter<int>("hltHadIsolDepth"); 
00087 
00088   calHLTHcalIsol_ = conf.getParameter<bool>("calHLTHcalIsol");
00089   calHLTEmIsol_ = conf.getParameter<bool>("calHLTEmIsol");
00090   calHLTEleTrkIsol_ = conf.getParameter<bool>("calHLTEleTrkIsol");
00091   calHLTPhoTrkIsol_ = conf.getParameter<bool>("calHLTPhoTrkIsol");
00092  
00093   trigCutParams_ = conf.getParameter<std::vector<edm::ParameterSet> >("triggerCuts"); //setupTriggers used to be in this function but had to be moved due to HLTConfigChanges (has to be called beginRun) so we have to save this for later.
00094 
00095   hltEleTrkIsolAlgo_ = new EgammaHLTTrackIsolation(hltEleTrkIsolPtMin_,hltEleTrkIsolOuterCone_,hltEleTrkIsolZSpan_,hltEleTrkIsolRSpan_,hltEleTrkIsolInnerCone_);
00096   hltPhoTrkIsolAlgo_ = new EgammaHLTTrackIsolation(hltPhoTrkIsolPtMin_,hltPhoTrkIsolOuterCone_,hltPhoTrkIsolZSpan_,hltPhoTrkIsolRSpan_,hltPhoTrkIsolInnerCone_);
00097                                                
00098 
00099 }
00100 
00101 //this code was taken out of OffHelper::setup due to HLTConfigProvider changes
00102 //it still assumes that this is called only once
00103 void OffHelper::setupTriggers(const HLTConfigProvider& hltConfig,const std::vector<std::string>& hltFiltersUsed)
00104 {
00105   hltFiltersUsed_ = hltFiltersUsed; //expensive but only do this once and faster ways could make things less clear
00106   //now work out how many objects are requires to pass filter for it to accept
00107   hltFiltersUsedWithNrCandsCut_.clear();
00108   for(size_t filterNr=0;filterNr<hltFiltersUsed_.size();filterNr++){
00109     hltFiltersUsedWithNrCandsCut_.push_back(std::make_pair(hltFiltersUsed_[filterNr],egHLT::trigTools::getMinNrObjsRequiredByFilter(hltFiltersUsed_[filterNr])));
00110   }
00111 
00112   //now loading the cuts for every trigger into our vector which stores them
00113   //only load cuts for triggers that are in hltFiltersUsed
00114   
00115   for(size_t trigNr=0;trigNr<trigCutParams_.size();trigNr++) {
00116     std::string trigName = trigCutParams_[trigNr].getParameter<std::string>("trigName");
00117     if(std::find(hltFiltersUsed_.begin(),hltFiltersUsed_.end(),trigName)!=hltFiltersUsed_.end()){ //perhaps I should sort hltFiltersUsed_....
00118       trigCuts_.push_back(std::make_pair(TrigCodes::getCode(trigName),OffEgSel(trigCutParams_[trigNr])));
00119       //   std::cout<<trigName<<std::endl<<"between"<<std::endl<<trigCutParams_[trigNr]<<std::endl<<"after"<<std::endl;
00120     }
00121   }
00122   trigCutParams_.clear();//dont need it any more, get rid of it
00123 
00124   //to make my life difficult, the scaled l1 paths are special
00125   //and arent stored in trigger event
00126   //to I have to figure out the path, see if it passes
00127   //and then hunt down the l1 seed filter and use that to match to the pho/ele
00128   //matching on l1 seed filter is not enough as that will be passed for normal 
00129   //electron triggers even if pre-scale hasnt fired
00130   l1PreScaledFilters_.clear();
00131   l1PreScaledPaths_.clear();
00132   l1PreAndSeedFilters_.clear();
00133   for(size_t filterNr=0;filterNr<hltFiltersUsed_.size();filterNr++){
00134     if(hltFiltersUsed_[filterNr].find("hltPreL1")==0){ //l1 prescaled path
00135       l1PreScaledFilters_.push_back(hltFiltersUsed_[filterNr]);
00136     }
00137   }
00138 
00139   egHLT::trigTools::translateFiltersToPathNames(hltConfig,l1PreScaledFilters_,l1PreScaledPaths_);
00140   if(l1PreScaledPaths_.size()==l1PreScaledFilters_.size()){
00141     for(size_t pathNr=0;pathNr<l1PreScaledPaths_.size();pathNr++){
00142      
00143       std::string l1SeedFilter =egHLT::trigTools::getL1SeedFilterOfPath(hltConfig,l1PreScaledPaths_[pathNr]);
00144       //---Morse====
00145       //std::cout<<l1PreScaledFilters_[pathNr]<<"  "<<l1PreScaledPaths_[pathNr]<<"  "<<l1SeedFilter<<std::endl;
00146       //------------
00147       l1PreAndSeedFilters_.push_back(std::make_pair(l1PreScaledFilters_[pathNr],l1SeedFilter));
00148     }
00149   }
00150 }
00151 
00152 int OffHelper::makeOffEvt(const edm::Event& edmEvent,const edm::EventSetup& setup,egHLT::OffEvt& offEvent)
00153 {
00154   offEvent.clear();
00155   int errCode=0; //excution stops as soon as an error is flagged
00156   if(errCode==0) errCode = getHandles(edmEvent,setup);
00157   if(errCode==0) errCode = fillOffEleVec(offEvent.eles());
00158   if(errCode==0) errCode = fillOffPhoVec(offEvent.phos());
00159   if(errCode==0) errCode = setTrigInfo(edmEvent, offEvent);
00160   if(errCode==0) offEvent.setJets(recoJets_);
00161   return errCode;
00162 }
00163 
00164 
00165 int OffHelper::getHandles(const edm::Event& event,const edm::EventSetup& setup)
00166 {
00167   try { 
00168     setup.get<CaloGeometryRecord>().get(caloGeom_);
00169     setup.get<CaloTopologyRecord>().get(caloTopology_);
00170     //setup.get<EcalSeverityLevelAlgoRcd>().get(ecalSeverityLevel_);
00171   }catch(...){
00172     return errCodes::Geom;
00173   }
00174   try {
00175     setup.get<IdealMagneticFieldRecord>().get(magField_);
00176   }catch(...){
00177     return errCodes::MagField;
00178   }
00179 
00180   //get objects
00181   if(!getHandle(event,triggerSummaryLabel_,trigEvt_)) return errCodes::TrigEvent; //must have this, otherwise skip event
00182   if(!getHandle(event,trigResultsTag_,trigResults_)) return errCodes::TrigEvent; //re using bit to minimise bug fix code changes
00183   if(!getHandle(event,electronsTag_,recoEles_)) return errCodes::OffEle; //need for electrons
00184   if(!getHandle(event,photonsTag_, recoPhos_)) return errCodes::OffPho; //need for photons
00185   if(!getHandle(event,caloJetsTag_,recoJets_)) return errCodes::OffJet; //need for electrons and photons
00186   if(!getHandle(event,vertexTag_,recoVertices_)) return errCodes::OffVertex; //need for eff vs nVertex
00187 
00188   //need for HLT isolations (rec hits also need for sigmaIPhiIPhi (ele/pho) and r9 pho)
00189   if(!getHandle(event,ecalRecHitsEBTag_,ebRecHits_)) return errCodes::EBRecHits;
00190   if(!getHandle(event,ecalRecHitsEETag_,eeRecHits_)) return errCodes::EERecHits;
00191   if(!getHandle(event,isolTrkTag_,isolTrks_)) return errCodes::IsolTrks;
00192   if(!getHandle(event,hbheHitsTag_, hbheHits_)) return errCodes::HBHERecHits; //I dont think we need hbhe rec-hits any more
00193   if(!getHandle(event,hfHitsTag_, hfHits_)) return errCodes::HFRecHits;//I dont think we need hf rec-hits any more
00194   if(!getHandle(event,beamSpotTag_,beamSpot_)) return errCodes::BeamSpot;
00195   if(!getHandle(event,caloTowersTag_,caloTowers_)) return errCodes::CaloTowers;
00196 
00197   
00198   return 0;
00199 }
00200 
00201 //this function coverts GsfElectrons to a format which is actually useful to me
00202 int OffHelper::fillOffEleVec(std::vector<OffEle>& egHLTOffEles)
00203 {
00204   egHLTOffEles.clear();
00205   egHLTOffEles.reserve(recoEles_->size());
00206   for(reco::GsfElectronCollection::const_iterator gsfIter=recoEles_->begin(); gsfIter!=recoEles_->end();++gsfIter){
00207     if(!gsfIter->ecalDrivenSeed()) continue; //avoid PF electrons (this is Eg HLT validation and HLT is ecal driven)
00208 
00209     int nVertex=0;
00210     for(reco::VertexCollection::const_iterator nVit=recoVertices_->begin(); nVit!=recoVertices_->end();++nVit){
00211       if( !nVit->isFake() 
00212           && nVit->ndof()>4  
00213           && std::fabs( nVit->z()<24.0) 
00214           && sqrt(nVit->x()*nVit->x() + nVit->y()*nVit->y())<2.0){nVertex++;}
00215     }
00216     //if(nVertex>20)std::cout<<"nVertex: "<<nVertex<<std::endl;
00217     OffEle::EventData eventData;
00218     eventData.NVertex=nVertex;
00219 
00220     OffEle::IsolData isolData;   
00221     fillIsolData(*gsfIter,isolData);
00222     
00223     OffEle::ClusShapeData clusShapeData;
00224     fillClusShapeData(*gsfIter,clusShapeData);
00225 
00226     OffEle::HLTData hltData;
00227     fillHLTData(*gsfIter,hltData);
00228 
00229     egHLTOffEles.push_back(OffEle(*gsfIter,clusShapeData,isolData,hltData,eventData));
00230     
00231     //now we would like to set the cut results
00232     OffEle& ele =  egHLTOffEles.back();
00233     ele.setCutCode(eleCuts_.getCutCode(ele));
00234     ele.setLooseCutCode(eleLooseCuts_.getCutCode(ele));
00235     
00236     std::vector<std::pair<TrigCodes::TrigBitSet,int> >trigCutsCutCodes;
00237     for(size_t i=0;i<trigCuts_.size();i++) trigCutsCutCodes.push_back(std::make_pair(trigCuts_[i].first,trigCuts_[i].second.getCutCode(ele)));
00238     ele.setTrigCutsCutCodes(trigCutsCutCodes);
00239   }//end loop over gsf electron collection
00240   return 0;
00241 }
00242 
00243 void OffHelper::fillIsolData(const reco::GsfElectron& ele,OffEle::IsolData& isolData)
00244 {
00245   EgammaTowerIsolation hcalIsolAlgo(hltHadIsolOuterCone_,hltHadIsolInnerCone_,hltHadIsolEtMin_,hltHadIsolDepth_,caloTowers_.product());
00246   EcalRecHitMetaCollection ebHits(*ebRecHits_);
00247   EcalRecHitMetaCollection eeHits(*eeRecHits_);
00248   EgammaRecHitIsolation ecalIsolAlgoEB(hltEMIsolOuterCone_,hltEMIsolInnerConeEB_,hltEMIsolEtaSliceEB_,
00249                                        hltEMIsolEtMinEB_,hltEMIsolEMinEB_,caloGeom_,&ebHits,ecalSeverityLevel_.product(),DetId::Ecal);
00250   EgammaRecHitIsolation ecalIsolAlgoEE(hltEMIsolOuterCone_,hltEMIsolInnerConeEE_,hltEMIsolEtaSliceEE_,
00251                                        hltEMIsolEtMinEE_,hltEMIsolEMinEE_,caloGeom_,&eeHits,ecalSeverityLevel_.product(),DetId::Ecal);
00252   
00253   isolData.ptTrks=ele.dr03TkSumPt();
00254   isolData.nrTrks=999; //no longer supported
00255   isolData.em= ele.dr03EcalRecHitSumEt();
00256   isolData.hadDepth1 =  ele.dr03HcalDepth1TowerSumEt();
00257   isolData.hadDepth2 =  ele.dr03HcalDepth2TowerSumEt();   
00258 
00259   //now time to do the HLT algos
00260   if(calHLTHcalIsol_) isolData.hltHad=hcalIsolAlgo.getTowerESum(&ele);
00261   else isolData.hltHad = 0.;
00262   if(calHLTEleTrkIsol_) isolData.hltTrksEle=hltEleTrkIsolAlgo_->electronPtSum(&(*(ele.gsfTrack())),isolTrks_.product());
00263   else isolData.hltTrksEle = 0.;
00264   if(calHLTPhoTrkIsol_){
00265     if(hltPhoTrkIsolCountTrks_) isolData.hltTrksPho=hltPhoTrkIsolAlgo_->photonTrackCount(&ele,isolTrks_.product(),false);
00266     else isolData.hltTrksPho=hltPhoTrkIsolAlgo_->photonPtSum(&ele,isolTrks_.product(),false);
00267   }
00268   else isolData.hltTrksPho = 0.;
00269   if(calHLTEmIsol_) isolData.hltEm = ecalIsolAlgoEB.getEtSum(&ele) + 
00270                                      ecalIsolAlgoEE.getEtSum(&ele);
00271   else isolData.hltEm = 0.;
00272   
00273 }
00274 
00275 
00276 void OffHelper::fillClusShapeData(const reco::GsfElectron& ele,OffEle::ClusShapeData& clusShapeData)
00277 {
00278   clusShapeData.sigmaEtaEta = ele.sigmaEtaEta();
00279   clusShapeData.sigmaIEtaIEta = ele.sigmaIetaIeta();
00280   double e5x5 = ele.e5x5();
00281   if(e5x5!=0.){
00282     clusShapeData.e1x5Over5x5 = ele.e1x5()/e5x5;
00283     clusShapeData.e2x5MaxOver5x5 = ele.e2x5Max()/e5x5;
00284   }else{
00285     clusShapeData.e1x5Over5x5 = -1;
00286     clusShapeData.e2x5MaxOver5x5 = -1;
00287   }
00288   
00289   //want to calculate r9, sigmaPhiPhi and sigmaIPhiIPhi, have to do old fashioned way
00290   const reco::BasicCluster& seedClus = *(ele.superCluster()->seed());
00291   const DetId seedDetId = seedClus.hitsAndFractions()[0].first; //note this may not actually be the seed hit but it doesnt matter because all hits will be in the barrel OR endcap
00292   if(seedDetId.subdetId()==EcalBarrel){
00293     std::vector<float> stdCov = EcalClusterTools::covariances(seedClus,ebRecHits_.product(),caloTopology_.product(),caloGeom_.product());
00294     std::vector<float> crysCov = EcalClusterTools::localCovariances(seedClus,ebRecHits_.product(),caloTopology_.product());
00295     clusShapeData.sigmaPhiPhi = sqrt(stdCov[2]);
00296       clusShapeData.sigmaIPhiIPhi =  sqrt(crysCov[2]);
00297       if(ele.superCluster()->rawEnergy()!=0.){
00298         clusShapeData.r9 = EcalClusterTools::e3x3(seedClus,ebRecHits_.product(),caloTopology_.product()) / ele.superCluster()->rawEnergy();
00299       }else clusShapeData.r9 = -1.;
00300        
00301   }else{
00302     std::vector<float> stdCov = EcalClusterTools::covariances(seedClus,eeRecHits_.product(),caloTopology_.product(),caloGeom_.product()); 
00303     std::vector<float> crysCov = EcalClusterTools::localCovariances(seedClus,eeRecHits_.product(),caloTopology_.product());
00304     clusShapeData.sigmaPhiPhi = sqrt(stdCov[2]);
00305     clusShapeData.sigmaIPhiIPhi = sqrt(crysCov[2]);
00306     if(ele.superCluster()->rawEnergy()!=0.){
00307         clusShapeData.r9 = EcalClusterTools::e3x3(seedClus,eeRecHits_.product(),caloTopology_.product()) / ele.superCluster()->rawEnergy();
00308     }else clusShapeData.r9 = -1.;
00309   } 
00310 }
00311 
00312 //reco approximations of hlt quantities
00313 void OffHelper::fillHLTData(const reco::GsfElectron& ele,OffEle::HLTData& hltData)
00314 {
00315   if(ele.closestCtfTrackRef().isNonnull() && 
00316      ele.closestCtfTrackRef()->extra().isNonnull()){
00317     reco::TrackRef ctfTrack = ele.closestCtfTrackRef();
00318     reco::SuperClusterRef scClus = ele.superCluster();
00319 
00320     //dEta
00321     const reco::BeamSpot::Point& bsPos = beamSpot_->position();     
00322     math::XYZPoint scPosWRTVtx(scClus->x()-bsPos.x(), scClus->y()-bsPos.y() , scClus->z()-ctfTrack->vz());
00323     hltData.dEtaIn = fabs(scPosWRTVtx.eta()-ctfTrack->eta());
00324 
00325     //dPhi: lifted straight from hlt code
00326     float deltaPhi=fabs(ctfTrack->outerPosition().phi()-scClus->phi());
00327     if(deltaPhi>6.283185308) deltaPhi -= 6.283185308;
00328     if(deltaPhi>3.141592654) deltaPhi = 6.283185308-deltaPhi;
00329     hltData.dPhiIn = deltaPhi;
00330     
00331     //invEInvP
00332     if(ele.ecalEnergy()!=0 && ctfTrack->p()!=0) hltData.invEInvP= 1/ele.ecalEnergy() - 1/ctfTrack->p();
00333     else hltData.invEInvP = 0;
00334   }else{
00335     hltData.dEtaIn =999;
00336     hltData.dPhiIn =999;
00337     hltData.invEInvP = 999;
00338 
00339   }
00340 }
00341 
00342 //this function coverts Photons to a format which more useful to me
00343 int OffHelper::fillOffPhoVec(std::vector<OffPho>& egHLTOffPhos)
00344 {
00345   egHLTOffPhos.clear();
00346   egHLTOffPhos.reserve(recoPhos_->size());
00347   for(reco::PhotonCollection::const_iterator phoIter=recoPhos_->begin(); phoIter!=recoPhos_->end();++phoIter){
00348 
00349     OffPho::IsolData isolData;  
00350     OffPho::ClusShapeData clusShapeData;
00351   
00352     fillIsolData(*phoIter,isolData);
00353     fillClusShapeData(*phoIter,clusShapeData);
00354  
00355     egHLTOffPhos.push_back(OffPho(*phoIter,clusShapeData,isolData));
00356     OffPho& pho =  egHLTOffPhos.back();
00357     pho.setCutCode(phoCuts_.getCutCode(pho));
00358     pho.setLooseCutCode(phoLooseCuts_.getCutCode(pho));
00359 
00360     std::vector<std::pair<TrigCodes::TrigBitSet,int> >trigCutsCutCodes;
00361     for(size_t i=0;i<trigCuts_.size();i++) trigCutsCutCodes.push_back(std::make_pair(trigCuts_[i].first,trigCuts_[i].second.getCutCode(pho)));
00362     pho.setTrigCutsCutCodes(trigCutsCutCodes); 
00363 
00364 
00365   }//end loop over photon collection
00366   return 0;
00367 }
00368 
00369 
00370 void OffHelper::fillIsolData(const reco::Photon& pho,OffPho::IsolData& isolData)
00371 {
00372   EgammaTowerIsolation hcalIsolAlgo(hltHadIsolOuterCone_,hltHadIsolInnerCone_,hltHadIsolEtMin_,hltHadIsolDepth_,caloTowers_.product());
00373   EcalRecHitMetaCollection ebHits(*ebRecHits_);
00374   EcalRecHitMetaCollection eeHits(*ebRecHits_);
00375   EgammaRecHitIsolation ecalIsolAlgoEB(hltEMIsolOuterCone_,hltEMIsolInnerConeEB_,hltEMIsolEtaSliceEB_,
00376                                        hltEMIsolEtMinEB_,hltEMIsolEMinEB_,caloGeom_,&ebHits,ecalSeverityLevel_.product(),DetId::Ecal);
00377   EgammaRecHitIsolation ecalIsolAlgoEE(hltEMIsolOuterCone_,hltEMIsolInnerConeEE_,hltEMIsolEtaSliceEE_,
00378                                        hltEMIsolEtMinEE_,hltEMIsolEMinEE_,caloGeom_,&eeHits,ecalSeverityLevel_.product(),DetId::Ecal);
00379   
00380   isolData.nrTrks = pho.nTrkHollowConeDR03();
00381   isolData.ptTrks = pho.trkSumPtHollowConeDR03();
00382   isolData.em = pho.ecalRecHitSumEtConeDR03();
00383   isolData.had = pho.hcalTowerSumEtConeDR03();  
00384 
00385   //now calculate hlt algos
00386   if(calHLTHcalIsol_) isolData.hltHad=hcalIsolAlgo.getTowerESum(&pho);
00387   else isolData.hltHad = 0.;
00388   if(calHLTPhoTrkIsol_){
00389     if(hltPhoTrkIsolCountTrks_) isolData.hltTrks=hltPhoTrkIsolAlgo_->photonTrackCount(&pho,isolTrks_.product(),false);
00390     else isolData.hltTrks=hltPhoTrkIsolAlgo_->photonPtSum(&pho,isolTrks_.product(),false);
00391   }
00392   else isolData.hltTrks = 0.;
00393   if(calHLTEmIsol_) isolData.hltEm = ecalIsolAlgoEB.getEtSum(&pho) + 
00394                                      ecalIsolAlgoEE.getEtSum(&pho);
00395   else isolData.hltEm = 0.;
00396   
00397 }
00398 
00399 void OffHelper::fillClusShapeData(const reco::Photon& pho,OffPho::ClusShapeData& clusShapeData)
00400 {
00401   clusShapeData.sigmaEtaEta = pho.sigmaEtaEta();
00402   clusShapeData.sigmaIEtaIEta =  pho.sigmaIetaIeta();
00403   double e5x5 =  pho.e5x5();
00404   if(e5x5!=0.){ //even though it is almost impossible for this to be 0., this code can never ever crash under any situation
00405     clusShapeData.e1x5Over5x5 =  pho.e1x5()/e5x5;
00406     clusShapeData.e2x5MaxOver5x5 = pho.e2x5()/e5x5;
00407   }else{
00408     clusShapeData.e1x5Over5x5 = -1;
00409     clusShapeData.e2x5MaxOver5x5 = -1;
00410   }
00411   clusShapeData.r9 = pho.r9();
00412   
00413   //sigmaPhiPhi and sigmaIPhiIPhi are not in object (and nor should they be) so have to get them old fashioned way
00414   //need to figure out if its in the barrel or endcap
00415   //get the first hit of the cluster and figure out if its barrel or endcap 
00416   const reco::BasicCluster& seedClus = *(pho.superCluster()->seed());
00417   const DetId seedDetId = seedClus.hitsAndFractions()[0].first; //note this may not actually be the seed hit but it doesnt matter because all hits will be in the barrel OR endcap (it is also incredably inefficient as it getHitsByDetId passes the vector by value not reference
00418   if(seedDetId.subdetId()==EcalBarrel){
00419     std::vector<float> stdCov = EcalClusterTools::covariances(seedClus,ebRecHits_.product(),caloTopology_.product(),caloGeom_.product());
00420     std::vector<float> crysCov = EcalClusterTools::localCovariances(seedClus,ebRecHits_.product(),caloTopology_.product());
00421     clusShapeData.sigmaPhiPhi = sqrt(stdCov[2]);
00422     clusShapeData.sigmaIPhiIPhi =  sqrt(crysCov[2]);
00423     }else{
00424     std::vector<float> stdCov = EcalClusterTools::covariances(seedClus,eeRecHits_.product(),caloTopology_.product(),caloGeom_.product()); 
00425     std::vector<float> crysCov = EcalClusterTools::localCovariances(seedClus,eeRecHits_.product(),caloTopology_.product());
00426     
00427     clusShapeData.sigmaPhiPhi = sqrt(stdCov[2]);
00428     clusShapeData.sigmaIPhiIPhi = sqrt(crysCov[2]); 
00429   }
00430 }  
00431 
00432 int OffHelper::setTrigInfo(const edm::Event & edmEvent, egHLT::OffEvt& offEvent)
00433 {
00434   TrigCodes::TrigBitSet evtTrigBits = trigTools::getFiltersPassed(hltFiltersUsedWithNrCandsCut_,trigEvt_.product(),hltTag_);
00435   //the l1 prescale paths dont have a filter with I can figure out if it passed or failed with so have to use TriggerResults
00436   if(l1PreScaledPaths_.size()==l1PreScaledFilters_.size()){ //check to ensure both vectors have same number of events incase of screw ups     
00437     const edm::TriggerNames & triggerNames = edmEvent.triggerNames(*trigResults_);
00438     for(size_t pathNr=0;pathNr<l1PreScaledPaths_.size();pathNr++){ //now we have to check the prescaled l1 trigger paths
00439       unsigned int pathIndex = triggerNames.triggerIndex(l1PreScaledPaths_[pathNr]);
00440       if(pathIndex<trigResults_->size() && trigResults_->accept(pathIndex)){
00441         evtTrigBits |=TrigCodes::getCode(l1PreScaledFilters_[pathNr]);
00442       }    
00443     }
00444   }
00445 
00446   offEvent.setEvtTrigBits(evtTrigBits);
00447 
00448   trigTools::setFiltersObjPasses(offEvent.eles(),hltFiltersUsed_,l1PreAndSeedFilters_,evtTrigBits,trigEvt_.product(),hltTag_);
00449   trigTools::setFiltersObjPasses(offEvent.phos(),hltFiltersUsed_,l1PreAndSeedFilters_,evtTrigBits,trigEvt_.product(),hltTag_); 
00450   return 0;
00451 }