CMS 3D CMS Logo

/data/refman/pasoursint/CMSSW_5_3_8_patch3/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   //Now get HLT p4 from triggerobject
00342   trigTools::fillHLTposition(ele,hltData,hltFiltersUsed_,trigEvt_.product(),hltTag_);
00343   //trigTools::fillHLTposition(phos(),hltFiltersUsed_,l1PreAndSeedFilters_,evtTrigBits,trigEvt_.product(),hltTag_); 
00344 }
00345 
00346 
00347 void OffHelper::fillHLTDataPho(const reco::Photon& pho,OffPho::HLTData& hltData)
00348 {
00349   //Now get HLT p4 from triggerobject
00350   trigTools::fillHLTposition(pho,hltData, hltFiltersUsed_,trigEvt_.product(),hltTag_);
00351   //trigTools::fillHLTposition(phos(),hltFiltersUsed_,l1PreAndSeedFilters_,evtTrigBits,trigEvt_.product(),hltTag_); 
00352 }
00353 
00354 
00355 
00356 
00357 //this function coverts Photons to a format which more useful to me
00358 int OffHelper::fillOffPhoVec(std::vector<OffPho>& egHLTOffPhos)
00359 {
00360   egHLTOffPhos.clear();
00361   egHLTOffPhos.reserve(recoPhos_->size());
00362   for(reco::PhotonCollection::const_iterator phoIter=recoPhos_->begin(); phoIter!=recoPhos_->end();++phoIter){
00363 
00364     OffPho::IsolData isolData;  
00365     OffPho::ClusShapeData clusShapeData;
00366   
00367     fillIsolData(*phoIter,isolData);
00368     fillClusShapeData(*phoIter,clusShapeData);
00369    
00370     OffPho::HLTData hltData;
00371     fillHLTDataPho(*phoIter,hltData); 
00372 
00373     egHLTOffPhos.push_back(OffPho(*phoIter,clusShapeData,isolData,hltData));
00374     OffPho& pho =  egHLTOffPhos.back();
00375     pho.setCutCode(phoCuts_.getCutCode(pho));
00376     pho.setLooseCutCode(phoLooseCuts_.getCutCode(pho));
00377 
00378     std::vector<std::pair<TrigCodes::TrigBitSet,int> >trigCutsCutCodes;
00379     for(size_t i=0;i<trigCuts_.size();i++) trigCutsCutCodes.push_back(std::make_pair(trigCuts_[i].first,trigCuts_[i].second.getCutCode(pho)));
00380     pho.setTrigCutsCutCodes(trigCutsCutCodes); 
00381 
00382 
00383   }//end loop over photon collection
00384   return 0;
00385 }
00386 
00387 
00388 void OffHelper::fillIsolData(const reco::Photon& pho,OffPho::IsolData& isolData)
00389 {
00390   EgammaTowerIsolation hcalIsolAlgo(hltHadIsolOuterCone_,hltHadIsolInnerCone_,hltHadIsolEtMin_,hltHadIsolDepth_,caloTowers_.product());
00391   EcalRecHitMetaCollection ebHits(*ebRecHits_);
00392   EcalRecHitMetaCollection eeHits(*ebRecHits_);
00393   EgammaRecHitIsolation ecalIsolAlgoEB(hltEMIsolOuterCone_,hltEMIsolInnerConeEB_,hltEMIsolEtaSliceEB_,
00394                                        hltEMIsolEtMinEB_,hltEMIsolEMinEB_,caloGeom_,&ebHits,ecalSeverityLevel_.product(),DetId::Ecal);
00395   EgammaRecHitIsolation ecalIsolAlgoEE(hltEMIsolOuterCone_,hltEMIsolInnerConeEE_,hltEMIsolEtaSliceEE_,
00396                                        hltEMIsolEtMinEE_,hltEMIsolEMinEE_,caloGeom_,&eeHits,ecalSeverityLevel_.product(),DetId::Ecal);
00397   
00398   isolData.nrTrks = pho.nTrkHollowConeDR03();
00399   isolData.ptTrks = pho.trkSumPtHollowConeDR03();
00400   isolData.em = pho.ecalRecHitSumEtConeDR03();
00401   isolData.had = pho.hcalTowerSumEtConeDR03();  
00402 
00403   //now calculate hlt algos
00404   if(calHLTHcalIsol_) isolData.hltHad=hcalIsolAlgo.getTowerESum(&pho);
00405   else isolData.hltHad = 0.;
00406   if(calHLTPhoTrkIsol_){
00407     if(hltPhoTrkIsolCountTrks_) isolData.hltTrks=hltPhoTrkIsolAlgo_->photonTrackCount(&pho,isolTrks_.product(),false);
00408     else isolData.hltTrks=hltPhoTrkIsolAlgo_->photonPtSum(&pho,isolTrks_.product(),false);
00409   }
00410   else isolData.hltTrks = 0.;
00411   if(calHLTEmIsol_) isolData.hltEm = ecalIsolAlgoEB.getEtSum(&pho) + 
00412                                      ecalIsolAlgoEE.getEtSum(&pho);
00413   else isolData.hltEm = 0.;
00414   
00415 }
00416 
00417 void OffHelper::fillClusShapeData(const reco::Photon& pho,OffPho::ClusShapeData& clusShapeData)
00418 {
00419   clusShapeData.sigmaEtaEta = pho.sigmaEtaEta();
00420   clusShapeData.sigmaIEtaIEta =  pho.sigmaIetaIeta();
00421   double e5x5 =  pho.e5x5();
00422   if(e5x5!=0.){ //even though it is almost impossible for this to be 0., this code can never ever crash under any situation
00423     clusShapeData.e1x5Over5x5 =  pho.e1x5()/e5x5;
00424     clusShapeData.e2x5MaxOver5x5 = pho.e2x5()/e5x5;
00425   }else{
00426     clusShapeData.e1x5Over5x5 = -1;
00427     clusShapeData.e2x5MaxOver5x5 = -1;
00428   }
00429   clusShapeData.r9 = pho.r9();
00430   
00431   //sigmaPhiPhi and sigmaIPhiIPhi are not in object (and nor should they be) so have to get them old fashioned way
00432   //need to figure out if its in the barrel or endcap
00433   //get the first hit of the cluster and figure out if its barrel or endcap 
00434   const reco::BasicCluster& seedClus = *(pho.superCluster()->seed());
00435   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
00436   if(seedDetId.subdetId()==EcalBarrel){
00437     std::vector<float> stdCov = EcalClusterTools::covariances(seedClus,ebRecHits_.product(),caloTopology_.product(),caloGeom_.product());
00438     std::vector<float> crysCov = EcalClusterTools::localCovariances(seedClus,ebRecHits_.product(),caloTopology_.product());
00439     clusShapeData.sigmaPhiPhi = sqrt(stdCov[2]);
00440     clusShapeData.sigmaIPhiIPhi =  sqrt(crysCov[2]);
00441     }else{
00442     std::vector<float> stdCov = EcalClusterTools::covariances(seedClus,eeRecHits_.product(),caloTopology_.product(),caloGeom_.product()); 
00443     std::vector<float> crysCov = EcalClusterTools::localCovariances(seedClus,eeRecHits_.product(),caloTopology_.product());
00444     
00445     clusShapeData.sigmaPhiPhi = sqrt(stdCov[2]);
00446     clusShapeData.sigmaIPhiIPhi = sqrt(crysCov[2]); 
00447   }
00448 }  
00449 
00450 int OffHelper::setTrigInfo(const edm::Event & edmEvent, egHLT::OffEvt& offEvent)
00451 {
00452   TrigCodes::TrigBitSet evtTrigBits = trigTools::getFiltersPassed(hltFiltersUsedWithNrCandsCut_,trigEvt_.product(),hltTag_);
00453   //the l1 prescale paths dont have a filter with I can figure out if it passed or failed with so have to use TriggerResults
00454   if(l1PreScaledPaths_.size()==l1PreScaledFilters_.size()){ //check to ensure both vectors have same number of events incase of screw ups     
00455     const edm::TriggerNames & triggerNames = edmEvent.triggerNames(*trigResults_);
00456     for(size_t pathNr=0;pathNr<l1PreScaledPaths_.size();pathNr++){ //now we have to check the prescaled l1 trigger paths
00457       unsigned int pathIndex = triggerNames.triggerIndex(l1PreScaledPaths_[pathNr]);
00458       if(pathIndex<trigResults_->size() && trigResults_->accept(pathIndex)){
00459         evtTrigBits |=TrigCodes::getCode(l1PreScaledFilters_[pathNr]);
00460       }    
00461     }
00462   }
00463 
00464   offEvent.setEvtTrigBits(evtTrigBits);
00465 
00466   trigTools::setFiltersObjPasses(offEvent.eles(),hltFiltersUsed_,l1PreAndSeedFilters_,evtTrigBits,trigEvt_.product(),hltTag_);
00467   trigTools::setFiltersObjPasses(offEvent.phos(),hltFiltersUsed_,l1PreAndSeedFilters_,evtTrigBits,trigEvt_.product(),hltTag_); 
00468   return 0;
00469 }