CMS 3D CMS Logo

/data/refman/pasoursint/CMSSW_4_2_9_HLT1_bphpatch4/src/HLTriggerOffline/Egamma/src/EmDQM.cc

Go to the documentation of this file.
00001 
00002 //                    Header file for this                                    //
00004 #include "HLTriggerOffline/Egamma/interface/EmDQM.h"
00005 
00007 //                    Collaborating Class Header                              //
00009 #include "FWCore/Framework/interface/Frameworkfwd.h"
00010 #include "FWCore/ParameterSet/interface/ParameterSet.h"
00011 #include "FWCore/Framework/interface/MakerMacros.h"
00012 #include "DataFormats/HLTReco/interface/TriggerEventWithRefs.h"
00013 #include "DataFormats/RecoCandidate/interface/RecoEcalCandidate.h"
00014 #include "DataFormats/EgammaCandidates/interface/Electron.h"
00015 //#include "DataFormats/EgammaCandidates/interface/Photon.h"
00016 #include "DataFormats/L1Trigger/interface/L1EmParticle.h"
00017 #include "DataFormats/L1Trigger/interface/L1EmParticleFwd.h"
00018 //#include "SimDataFormats/HepMCProduct/interface/HepMCProduct.h"
00019 #include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h"
00020 #include "FWCore/MessageLogger/interface/MessageLogger.h"
00021 #include "DataFormats/Common/interface/AssociationMap.h"
00022 
00023 #include "DataFormats/Common/interface/Handle.h"
00024 #include "DataFormats/Common/interface/RefToBase.h"
00025 #include "FWCore/ServiceRegistry/interface/Service.h"
00026 #include "FWCore/Utilities/interface/Exception.h"
00027 #include "DataFormats/HLTReco/interface/TriggerTypeDefs.h"
00029 //                           Root include files                               //
00031 #include "TFile.h"
00032 #include "TDirectory.h"
00033 #include "TH1F.h"
00034 #include <iostream>
00035 #include <string>
00036 #include <Math/VectorUtil.h>
00037 using namespace ROOT::Math::VectorUtil ;
00038 
00039 
00041 //                             Constructor                                    //
00043 EmDQM::EmDQM(const edm::ParameterSet& pset)  
00044 {
00045 
00046   dbe = edm::Service < DQMStore > ().operator->();
00047   dbe->setVerbose(0);
00048 
00049 
00051   //          Read from configuration file                  //
00053   dirname_="HLT/HLTEgammaValidation/"+pset.getParameter<std::string>("@module_label");
00054   dbe->setCurrentFolder(dirname_);
00055 
00056   triggerobjwithrefs = pset.getParameter<edm::InputTag>("triggerobject");
00057   // paramters for generator study 
00058   reqNum    = pset.getParameter<unsigned int>("reqNum");
00059   pdgGen    = pset.getParameter<int>("pdgGen");
00060   genEtaAcc = pset.getParameter<double>("genEtaAcc");
00061   genEtAcc  = pset.getParameter<double>("genEtAcc");
00062   // plotting paramters (untracked because they don't affect the physics)
00063   plotPtMin  = pset.getUntrackedParameter<double>("PtMin",0.);
00064   plotPtMax  = pset.getUntrackedParameter<double>("PtMax",1000.);
00065   plotEtaMax = pset.getUntrackedParameter<double>("EtaMax", 2.7);
00066   plotPhiMax = pset.getUntrackedParameter<double>("PhiMax", 3.15);
00067   plotBins   = pset.getUntrackedParameter<unsigned int>("Nbins",40);
00068   plotMinEtForEtaEffPlot = pset.getUntrackedParameter<unsigned int>("minEtForEtaEffPlot", 15);
00069   useHumanReadableHistTitles = pset.getUntrackedParameter<bool>("useHumanReadableHistTitles", false);
00070 
00071   //preselction cuts 
00072   gencutCollection_= pset.getParameter<edm::InputTag>("cutcollection");
00073   gencut_          = pset.getParameter<int>("cutnum");
00074 
00076   //         Read in the Vector of Parameter Sets.          //
00077   //           Information for each filter-step             //
00079   std::vector<edm::ParameterSet> filters = 
00080        pset.getParameter<std::vector<edm::ParameterSet> >("filters");
00081 
00082   int i = 0;
00083   for(std::vector<edm::ParameterSet>::iterator filterconf = filters.begin() ; filterconf != filters.end() ; filterconf++)
00084   {
00085 
00086     theHLTCollectionLabels.push_back(filterconf->getParameter<edm::InputTag>("HLTCollectionLabels"));
00087     theHLTOutputTypes.push_back(filterconf->getParameter<int>("theHLTOutputTypes"));
00088     // Grab the human-readable name, if it is not specified, use the Collection Label
00089     theHLTCollectionHumanNames.push_back(filterconf->getUntrackedParameter<std::string>("HLTCollectionHumanName",theHLTCollectionLabels[i].label()));
00090 
00091     std::vector<double> bounds = filterconf->getParameter<std::vector<double> >("PlotBounds");
00092     // If the size of plot "bounds" vector != 2, abort
00093     assert(bounds.size() == 2);
00094     plotBounds.push_back(std::pair<double,double>(bounds[0],bounds[1]));
00095     isoNames.push_back(filterconf->getParameter<std::vector<edm::InputTag> >("IsoCollections"));
00096     // If the size of the isoNames vector is not greater than zero, abort
00097     assert(isoNames.back().size()>0);
00098     if (isoNames.back().at(0).label()=="none") {
00099       plotiso.push_back(false);
00100     } else {
00101       plotiso.push_back(true);
00102     }
00103     i++;
00104   } // END of loop over parameter sets
00105 
00106   // Record number of HLTCollectionLabels
00107   numOfHLTCollectionLabels = theHLTCollectionLabels.size();
00108   
00109 }
00110 
00111 
00113 //       method called once each job just before starting event loop          //
00115 void 
00116 EmDQM::beginJob()
00117 {
00118   //edm::Service<TFileService> fs;
00119   dbe->setCurrentFolder(dirname_);
00120 
00122   //  Set up Histogram of Effiency vs Step.                 //
00123   //   theHLTCollectionLabels is a vector of InputTags      //
00124   //    from the configuration file.                        //
00126 
00127   std::string histName="total_eff";
00128   std::string histTitle = "total events passing";
00129   // This plot will have bins equal to 2+(number of
00130   //        HLTCollectionLabels in the config file)
00131   total = dbe->book1D(histName.c_str(),histTitle.c_str(),numOfHLTCollectionLabels+2,0,numOfHLTCollectionLabels+2);
00132   total->setBinLabel(numOfHLTCollectionLabels+1,"Total");
00133   total->setBinLabel(numOfHLTCollectionLabels+2,"Gen");
00134   for (unsigned int u=0; u<numOfHLTCollectionLabels; u++){total->setBinLabel(u+1,theHLTCollectionLabels[u].label().c_str());}
00135 
00136   histName="total_eff_MC_matched";
00137   histTitle="total events passing (mc matched)";
00138   totalmatch = dbe->book1D(histName.c_str(),histTitle.c_str(),numOfHLTCollectionLabels+2,0,numOfHLTCollectionLabels+2);
00139   totalmatch->setBinLabel(numOfHLTCollectionLabels+1,"Total");
00140   totalmatch->setBinLabel(numOfHLTCollectionLabels+2,"Gen");
00141   for (unsigned int u=0; u<numOfHLTCollectionLabels; u++){totalmatch->setBinLabel(u+1,theHLTCollectionLabels[u].label().c_str());}
00142 
00143   MonitorElement* tmphisto;
00144   MonitorElement* tmpiso;
00145 
00147   // Set up generator-level histograms                      //
00149   std::string pdgIdString;
00150   switch(pdgGen) {
00151   case 11:
00152     pdgIdString="Electron";break;
00153   case 22:
00154     pdgIdString="Photon";break;
00155   default:
00156     pdgIdString="Particle";
00157   }
00158 
00159   histName = "gen_et";
00160   histTitle= "E_{T} of " + pdgIdString + "s" ;
00161   etgen =  dbe->book1D(histName.c_str(),histTitle.c_str(),plotBins,plotPtMin,plotPtMax);
00162   histName = "gen_eta";
00163   histTitle= "#eta of "+ pdgIdString +"s " ;
00164   etagen = dbe->book1D(histName.c_str(),histTitle.c_str(),plotBins,-plotEtaMax,plotEtaMax);
00165   histName = "gen_phi";
00166   histTitle= "#phi of "+ pdgIdString +"s " ;
00167   phigen = dbe->book1D(histName.c_str(),histTitle.c_str(),plotBins,-plotPhiMax,plotPhiMax);
00168 
00169   
00170 
00172   //  Set up histograms of HLT objects                      //
00174 
00175   // Determine what strings to use for histogram titles
00176   std::vector<std::string> HltHistTitle;
00177   if ( theHLTCollectionHumanNames.size() == numOfHLTCollectionLabels && useHumanReadableHistTitles ) {
00178     HltHistTitle = theHLTCollectionHumanNames;
00179   } else {
00180     for (unsigned int i =0; i < numOfHLTCollectionLabels; i++) {
00181       HltHistTitle.push_back(theHLTCollectionLabels[i].label());
00182     }
00183   }
00184  
00185   for(unsigned int i = 0; i< numOfHLTCollectionLabels ; i++){
00186     // Et distribution of HLT objects passing filter i
00187     histName = theHLTCollectionLabels[i].label()+"et_all";
00188     histTitle = HltHistTitle[i]+" Et (ALL)";
00189     tmphisto =  dbe->book1D(histName.c_str(),histTitle.c_str(),plotBins,plotPtMin,plotPtMax);
00190     ethist.push_back(tmphisto);
00191     
00192     // Eta distribution of HLT objects passing filter i
00193     histName = theHLTCollectionLabels[i].label()+"eta_all";
00194     histTitle = HltHistTitle[i]+" #eta (ALL)";
00195     tmphisto =  dbe->book1D(histName.c_str(),histTitle.c_str(),plotBins,-plotEtaMax,plotEtaMax);
00196     etahist.push_back(tmphisto);          
00197 
00198     // Phi distribution of HLT objects passing filter i
00199     histName = theHLTCollectionLabels[i].label()+"phi_all";
00200     histTitle = HltHistTitle[i]+" #phi (ALL)";
00201     tmphisto =  dbe->book1D(histName.c_str(),histTitle.c_str(),plotBins,-plotPhiMax,plotPhiMax);
00202     phihist.push_back(tmphisto);
00203 
00204 
00205     // Et distribution of gen object matching HLT object passing filter i
00206     histName = theHLTCollectionLabels[i].label()+"et_MC_matched";
00207     histTitle = HltHistTitle[i]+" Et (MC matched)";
00208     tmphisto =  dbe->book1D(histName.c_str(),histTitle.c_str(),plotBins,plotPtMin,plotPtMax);
00209     ethistmatch.push_back(tmphisto);
00210     
00211     // Eta distribution of gen object matching HLT object passing filter i
00212     histName = theHLTCollectionLabels[i].label()+"eta_MC_matched";
00213     histTitle = HltHistTitle[i]+" #eta (MC matched)";
00214     tmphisto =  dbe->book1D(histName.c_str(),histTitle.c_str(),plotBins,-plotEtaMax,plotEtaMax);
00215     etahistmatch.push_back(tmphisto);
00216 
00217     // Phi distribution of gen object matching HLT object passing filter i
00218     histName = theHLTCollectionLabels[i].label()+"phi_MC_matched";
00219     histTitle = HltHistTitle[i]+" #phi (MC matched)";
00220     tmphisto =  dbe->book1D(histName.c_str(),histTitle.c_str(),plotBins,-plotPhiMax,plotPhiMax);
00221     phihistmatch.push_back(tmphisto);
00222 
00223 
00224 
00225     // Et distribution of HLT object that is closest delta-R match to sorted gen particle(s)
00226     histName  = theHLTCollectionLabels[i].label()+"et";
00227     histTitle = HltHistTitle[i]+" Et";
00228     tmphisto  = dbe->book1D(histName.c_str(),histTitle.c_str(),plotBins,plotPtMin,plotPtMax);
00229     histEtOfHltObjMatchToGen.push_back(tmphisto);
00230 
00231     // eta distribution of HLT object that is closest delta-R match to sorted gen particle(s)
00232     histName  = theHLTCollectionLabels[i].label()+"eta";
00233     histTitle = HltHistTitle[i]+" eta";
00234     tmphisto  = dbe->book1D(histName.c_str(),histTitle.c_str(),plotBins,-plotEtaMax,plotEtaMax);
00235     histEtaOfHltObjMatchToGen.push_back(tmphisto);
00236 
00237     // phi distribution of HLT object that is closest delta-R match to sorted gen particle(s)
00238     histName  = theHLTCollectionLabels[i].label()+"phi";
00239     histTitle = HltHistTitle[i]+" phi";
00240     tmphisto  = dbe->book1D(histName.c_str(),histTitle.c_str(),plotBins,-plotPhiMax,plotPhiMax);
00241     histPhiOfHltObjMatchToGen.push_back(tmphisto);
00242 
00243 
00244 
00245     if (!plotiso[i]) {
00246       tmpiso = NULL;
00247       etahistiso.push_back(tmpiso);
00248       phihistiso.push_back(tmpiso);
00249       ethistiso.push_back(tmpiso);
00250       etahistisomatch.push_back(tmpiso);
00251       phihistisomatch.push_back(tmpiso);
00252       ethistisomatch.push_back(tmpiso);
00253       histEtaIsoOfHltObjMatchToGen.push_back(tmpiso);
00254       histPhiIsoOfHltObjMatchToGen.push_back(tmpiso);
00255       histEtIsoOfHltObjMatchToGen.push_back(tmpiso);
00256     } else {
00257       // 2D plot: Isolation values vs eta for all objects
00258       histName  = theHLTCollectionLabels[i].label()+"eta_isolation_all";
00259       histTitle = HltHistTitle[i]+" isolation vs #eta (all)";
00260       tmpiso    = dbe->book2D(histName.c_str(),histTitle.c_str(),plotBins,-plotEtaMax,plotEtaMax,plotBins,plotBounds[i].first,plotBounds[i].second);
00261       etahistiso.push_back(tmpiso);
00262 
00263       // 2D plot: Isolation values vs phi for all objects
00264       histName  = theHLTCollectionLabels[i].label()+"phi_isolation_all";
00265       histTitle = HltHistTitle[i]+" isolation vs #phi (all)";
00266       tmpiso    = dbe->book2D(histName.c_str(),histTitle.c_str(),plotBins,-plotPhiMax,plotPhiMax,plotBins,plotBounds[i].first,plotBounds[i].second);
00267       phihistiso.push_back(tmpiso);
00268 
00269 
00270       // 2D plot: Isolation values vs et for all objects
00271       histName  = theHLTCollectionLabels[i].label()+"et_isolation_all";
00272       histTitle = HltHistTitle[i]+" isolation vs Et (all)";
00273       tmpiso    = dbe->book2D(histName.c_str(),histTitle.c_str(),plotBins,plotPtMin,plotPtMax,plotBins,plotBounds[i].first,plotBounds[i].second);
00274       ethistiso.push_back(tmpiso);
00275 
00276       // 2D plot: Isolation values vs eta for matched objects
00277       histName  = theHLTCollectionLabels[i].label()+"eta_isolation_MC_matched";
00278       histTitle = HltHistTitle[i]+" isolation vs #eta (mc matched)";
00279       tmpiso    = dbe->book2D(histName.c_str(),histTitle.c_str(),plotBins,-plotEtaMax,plotEtaMax,plotBins,plotBounds[i].first,plotBounds[i].second);
00280       etahistisomatch.push_back(tmpiso);
00281 
00282       // 2D plot: Isolation values vs phi for matched objects
00283       histName  = theHLTCollectionLabels[i].label()+"phi_isolation_MC_matched";
00284       histTitle = HltHistTitle[i]+" isolation vs #phi (mc matched)";
00285       tmpiso    = dbe->book2D(histName.c_str(),histTitle.c_str(),plotBins,-plotPhiMax,plotPhiMax,plotBins,plotBounds[i].first,plotBounds[i].second);
00286       phihistisomatch.push_back(tmpiso);
00287 
00288 
00289       // 2D plot: Isolation values vs et for matched objects
00290       histName  = theHLTCollectionLabels[i].label()+"et_isolation_MC_matched";
00291       histTitle = HltHistTitle[i]+" isolation vs Et (mc matched)";
00292       tmpiso    = dbe->book2D(histName.c_str(),histTitle.c_str(),plotBins,plotPtMin,plotPtMax,plotBins,plotBounds[i].first,plotBounds[i].second);
00293       ethistisomatch.push_back(tmpiso);
00294 
00295 
00296       // 2D plot: Isolation values vs eta for HLT object that 
00297       // is closest delta-R match to sorted gen particle(s)
00298       histName  = theHLTCollectionLabels[i].label()+"eta_isolation";
00299       histTitle = HltHistTitle[i]+" isolation vs #eta";
00300       tmpiso    = dbe->book2D(histName.c_str(),histTitle.c_str(),plotBins,-plotEtaMax,plotEtaMax,plotBins,plotBounds[i].first,plotBounds[i].second);
00301       histEtaIsoOfHltObjMatchToGen.push_back(tmpiso);
00302 
00303       // 2D plot: Isolation values vs phi for HLT object that
00304       // is closest delta-R match to sorted gen particle(s)
00305       histName  = theHLTCollectionLabels[i].label()+"phi_isolation";
00306       histTitle = HltHistTitle[i]+" isolation vs #phi";
00307       tmpiso    = dbe->book2D(histName.c_str(),histTitle.c_str(),plotBins,-plotPhiMax,plotPhiMax,plotBins,plotBounds[i].first,plotBounds[i].second);
00308       histPhiIsoOfHltObjMatchToGen.push_back(tmpiso);
00309 
00310       // 2D plot: Isolation values vs et for HLT object that 
00311       // is closest delta-R match to sorted gen particle(s)
00312       histName  = theHLTCollectionLabels[i].label()+"et_isolation";
00313       histTitle = HltHistTitle[i]+" isolation vs Et";
00314       tmpiso    = dbe->book2D(histName.c_str(),histTitle.c_str(),plotBins,plotPtMin,plotPtMax,plotBins,plotBounds[i].first,plotBounds[i].second);
00315       histEtIsoOfHltObjMatchToGen.push_back(tmpiso);
00316     } // END of HLT histograms
00317 
00318   }
00319 }
00320 
00321 
00323 //                                Destructor                                  //
00325 EmDQM::~EmDQM(){
00326 }
00327 
00328 
00330 //                     method called to for each event                        //
00332 void 
00333 EmDQM::analyze(const edm::Event & event , const edm::EventSetup& setup)
00334 {
00335   
00337   //           Check if there's enough gen particles        //
00338   //             of interest                                //
00340   edm::Handle< edm::View<reco::Candidate> > cutCounter;
00341   event.getByLabel(gencutCollection_,cutCounter);
00342   if (cutCounter->size() < (unsigned int)gencut_) {
00343     //edm::LogWarning("EmDQM") << "Less than "<< gencut_ <<" gen particles with pdgId=" << pdgGen;
00344     return;
00345   }
00346 
00347 
00348   // fill L1 and HLT info
00349   // get objects possed by each filter
00350   edm::Handle<trigger::TriggerEventWithRefs> triggerObj;
00351   event.getByLabel(triggerobjwithrefs,triggerObj); 
00352   if(!triggerObj.isValid()) { 
00353     edm::LogWarning("EmDQM") << "triggerobjwithrefs=" << triggerobjwithrefs;
00354     return;
00355   }
00356 
00357 
00359   // Decide if this was an event of interest.               //
00360   //  Did the highest energy particles happen               //
00361   //  to have |eta| < 2.5 ?  Then continue.                 //
00363   edm::Handle< edm::View<reco::GenParticle> > genParticles;
00364   event.getByLabel("genParticles", genParticles);
00365   if(!genParticles.isValid()) { 
00366     edm::LogWarning("EmDQM") << "genParticles invalid.";
00367     return;
00368   }
00369 
00370   std::vector<reco::GenParticle> allSortedGenParticles;
00371 
00372   for(edm::View<reco::GenParticle>::const_iterator currentGenParticle = genParticles->begin(); currentGenParticle != genParticles->end(); currentGenParticle++){
00373 
00374     if (  !( abs((*currentGenParticle).pdgId())==pdgGen  && (*currentGenParticle).status()==1 && (*currentGenParticle).et() > 2.0)  )  continue;
00375 
00376     reco::GenParticle tmpcand( *(currentGenParticle) );
00377     allSortedGenParticles.push_back(tmpcand);
00378   }
00379 
00380   std::sort(allSortedGenParticles.begin(), allSortedGenParticles.end(),pTGenComparator_);
00381 
00382   // Were enough high energy gen particles found?
00383   if (allSortedGenParticles.size() < gencut_) {
00384     // if no, throw event away
00385     return;
00386   }
00387 
00388   // We now have a sorted collection of all generated particles
00389   // with pdgId = pdgGen.
00390   // Loop over them to see if the top gen particles have eta within acceptance
00391   bool keepEvent = true;
00392   for (unsigned int i = 0 ; i < gencut_ ; i++ ) {
00393     bool inECALgap = fabs(allSortedGenParticles[i].eta()) > 1.4442 && fabs(allSortedGenParticles[i].eta()) < 1.556;
00394     if ( (fabs(allSortedGenParticles[i].eta()) > genEtaAcc) || inECALgap ) {
00395       //edm::LogWarning("EmDQM") << "Throwing event away. Gen particle with pdgId="<< allSortedGenParticles[i].pdgId() <<"; et="<< allSortedGenParticles[i].et() <<"; and eta="<< allSortedGenParticles[i].eta() <<" beyond acceptance.";
00396       keepEvent=false;
00397       break;
00398     }
00399   }
00400   if (!keepEvent) {
00401     return;
00402   }
00403 
00404 
00405   // It was an event worth keeping. Continue.
00406 
00408   //  Fill the bin labeled "Total"                          //
00409   //   This will be the number of events looked at.         //
00411   total->Fill(numOfHLTCollectionLabels+0.5);
00412   totalmatch->Fill(numOfHLTCollectionLabels+0.5);
00413 
00414 
00416   //               Fill generator info                      //
00418   // the gencut_ highest Et generator objects of the preselected type are our matches
00419 
00420   std::vector<reco::Particle> sortedGen;
00421   for(edm::View<reco::Candidate>::const_iterator genpart = cutCounter->begin(); genpart != cutCounter->end();genpart++){
00422     reco::Particle tmpcand(  genpart->charge(), genpart->p4(), genpart->vertex(),genpart->pdgId(),genpart->status() );
00423     sortedGen.push_back(tmpcand);
00424   }
00425   std::sort(sortedGen.begin(),sortedGen.end(),pTComparator_ );
00426 
00427   // Now the collection of gen particles is sorted by pt.
00428   // So, remove all particles from the collection so that we 
00429   // only have the top "1 thru gencut_" particles in it
00430   sortedGen.erase(sortedGen.begin()+gencut_,sortedGen.end());
00431   if (gencut_ != sortedGen.size() ){
00432     return;
00433   }
00434 
00435 
00436   for (unsigned int i = 0 ; i < gencut_ ; i++ ) {
00437     etgen ->Fill( sortedGen[i].et()  ); //validity has been implicitily checked by the cut on gencut_ above
00438     if (sortedGen[i].et() > plotMinEtForEtaEffPlot) {
00439       etagen->Fill( sortedGen[i].eta() );
00440       phigen->Fill( sortedGen[i].phi() );
00441     }
00442   } // END of loop over Generated particles
00443   if (gencut_ >= reqNum) total->Fill(numOfHLTCollectionLabels+1.5); // this isn't really needed anymore keep for backward comp.
00444   if (gencut_ >= reqNum) totalmatch->Fill(numOfHLTCollectionLabels+1.5); // this isn't really needed anymore keep for backward comp.
00445           
00446 
00448   //            Loop over filter modules                    //
00450   for(unsigned int n=0; n < numOfHLTCollectionLabels ; n++) {
00451     // These numbers are from the Parameter Set, such as:
00452     //   theHLTOutputTypes = cms.uint32(100)
00453     switch(theHLTOutputTypes[n]) 
00454     {
00455       case trigger::TriggerL1NoIsoEG: // Non-isolated Level 1
00456         fillHistos<l1extra::L1EmParticleCollection>(triggerObj,event,n,sortedGen);break;
00457       case trigger::TriggerL1IsoEG: // Isolated Level 1
00458         fillHistos<l1extra::L1EmParticleCollection>(triggerObj,event,n,sortedGen);break;
00459       case trigger::TriggerPhoton: // Photon 
00460         fillHistos<reco::RecoEcalCandidateCollection>(triggerObj,event,n,sortedGen);break;
00461       case trigger::TriggerElectron: // Electron 
00462         fillHistos<reco::ElectronCollection>(triggerObj,event,n,sortedGen);break;
00463       case trigger::TriggerCluster: // TriggerCluster
00464         fillHistos<reco::RecoEcalCandidateCollection>(triggerObj,event,n,sortedGen);break;
00465       default: 
00466         throw(cms::Exception("Release Validation Error") << "HLT output type not implemented: theHLTOutputTypes[n]" );
00467     }
00468   } // END of loop over filter modules
00469 }
00470 
00471 
00473 // fillHistos                                                                 //
00474 //   Called by analyze method.                                                //
00476 template <class T> void EmDQM::fillHistos(edm::Handle<trigger::TriggerEventWithRefs>& triggerObj,const edm::Event& iEvent ,unsigned int n,std::vector<reco::Particle>& sortedGen)
00477 {
00478   std::vector<edm::Ref<T> > recoecalcands;
00479   if ( ( triggerObj->filterIndex(theHLTCollectionLabels[n])>=triggerObj->size() )){ // only process if available
00480     return;
00481   }
00482 
00484   //      Retrieve saved filter objects                     //
00486   triggerObj->getObjects(triggerObj->filterIndex(theHLTCollectionLabels[n]),theHLTOutputTypes[n],recoecalcands);
00487   //Danger: special case, L1 non-isolated
00488   // needs to be merged with L1 iso
00489   if (theHLTOutputTypes[n] == trigger::TriggerL1NoIsoEG){
00490     std::vector<edm::Ref<T> > isocands;
00491     triggerObj->getObjects(triggerObj->filterIndex(theHLTCollectionLabels[n]),trigger::TriggerL1IsoEG,isocands);
00492     if (isocands.size()>0) 
00493       {
00494         for (unsigned int i=0; i < isocands.size(); i++)
00495           recoecalcands.push_back(isocands[i]);
00496       }
00497   } // END of if theHLTOutputTypes == 82
00498   
00499 
00500   if (recoecalcands.size() < 1){ // stop if no object passed the previous filter
00501     return;
00502   }
00503 
00504 
00505   if (recoecalcands.size() >= reqNum ) 
00506     total->Fill(n+0.5);
00507 
00508 
00510   // check for validity                            //
00511   // prevents crash in CMSSW_3_1_0_pre6            //
00513   for (unsigned int j=0; j<recoecalcands.size(); j++){
00514     if(!( recoecalcands.at(j).isAvailable())){
00515       edm::LogError("EmDQM") << "Event content inconsistent: TriggerEventWithRefs contains invalid Refs" << std::endl << "invalid refs for: " << theHLTCollectionLabels[n].label();
00516       return;
00517     }
00518   }
00519 
00521   // Loop over the Generated Particles, and find the        //
00522   // closest HLT object match.                              //
00524   for (unsigned int i=0; i < gencut_; i++) {
00525     math::XYZVector currentGenParticleMomentum = sortedGen[i].momentum();
00526 
00527     float closestDeltaR = 0.5;
00528     int closestEcalCandIndex = -1;
00529     for (unsigned int j=0; j<recoecalcands.size(); j++) {
00530       float deltaR = DeltaR(recoecalcands[j]->momentum(),currentGenParticleMomentum);
00531 
00532       if (deltaR < closestDeltaR) {
00533         closestDeltaR = deltaR;
00534         closestEcalCandIndex = j;
00535       }
00536     }
00537 
00538     // If an HLT object was found within some delta-R
00539     // of this gen particle, store it in a histogram
00540     if ( closestEcalCandIndex >= 0 ) {
00541       histEtOfHltObjMatchToGen[n] ->Fill( recoecalcands[closestEcalCandIndex]->et()  );
00542       histEtaOfHltObjMatchToGen[n]->Fill( recoecalcands[closestEcalCandIndex]->eta() );
00543       histPhiOfHltObjMatchToGen[n]->Fill( recoecalcands[closestEcalCandIndex]->phi() );
00544       
00545       // Also store isolation info
00546       if (n+1 < numOfHLTCollectionLabels){ // can't plot beyond last
00547         if (plotiso[n+1] ){  // only plot if requested in config
00548           for (unsigned int j =  0 ; j < isoNames[n+1].size() ;j++  ){
00549             edm::Handle<edm::AssociationMap<edm::OneToValue< T , float > > > depMap; 
00550             iEvent.getByLabel(isoNames[n+1].at(j),depMap);
00551             if (depMap.isValid()){ //Map may not exist if only one candidate passes a double filter
00552               typename edm::AssociationMap<edm::OneToValue< T , float > >::const_iterator mapi = depMap->find(recoecalcands[closestEcalCandIndex]);
00553               if (mapi!=depMap->end()) {  // found candidate in isolation map! 
00554                 histEtaIsoOfHltObjMatchToGen[n+1]->Fill( recoecalcands[closestEcalCandIndex]->eta(),mapi->val);
00555                 histPhiIsoOfHltObjMatchToGen[n+1]->Fill( recoecalcands[closestEcalCandIndex]->phi(),mapi->val);
00556                 histEtIsoOfHltObjMatchToGen[n+1] ->Fill( recoecalcands[closestEcalCandIndex]->et(), mapi->val);
00557               }
00558             }
00559           }
00560         }
00561       }
00562     } // END of if closestEcalCandIndex >= 0
00563   }
00564 
00565 
00566 
00568   //  Loop over all HLT objects in this filter step, and    //
00569   //  fill histograms.                                      //
00571   //  bool foundAllMatches = false;
00572   //  unsigned int numOfHLTobjectsMatched = 0;
00573   for (unsigned int i=0; i<recoecalcands.size(); i++) {
00574     /*
00575     // See if this HLT object has a gen-level match
00576     float closestGenParticleDr = 99.0;
00577     for(unsigned int j =0; j < gencut_; j++) {
00578       math::XYZVector currentGenParticle = sortedGen[j].momentum();
00579 
00580       double currentDeltaR = DeltaR(recoecalcands[i]->momentum(),currentGenParticle);
00581       if ( currentDeltaR < closestGenParticleDr ) {
00582         closestGenParticleDr = currentDeltaR;
00583       }
00584     }
00585     // If this HLT object did not have a gen particle match, go to next HLT object
00586     if ( !(fabs(closestGenParticleDr < 0.3)) ) continue;
00587  
00588     numOfHLTobjectsMatched++;
00589     if (numOfHLTobjectsMatched >= gencut_) foundAllMatches=true;
00590     */
00591     // Fill HLT object histograms
00592     ethist[n] ->Fill(recoecalcands[i]->et() );
00593     etahist[n]->Fill(recoecalcands[i]->eta() );
00594     phihist[n]->Fill(recoecalcands[i]->phi() );
00595 
00597     //  Plot isolation variables (show the not-yet-cut        //
00598     //  isolation, i.e. associated to next filter)            //
00600     if ( n+1 < numOfHLTCollectionLabels ) { // can't plot beyond last
00601       if (plotiso[n+1]) {
00602         for (unsigned int j =  0 ; j < isoNames[n+1].size() ;j++  ){
00603           edm::Handle<edm::AssociationMap<edm::OneToValue< T , float > > > depMap; 
00604           iEvent.getByLabel(isoNames[n+1].at(j),depMap);
00605           if (depMap.isValid()){ //Map may not exist if only one candidate passes a double filter
00606             typename edm::AssociationMap<edm::OneToValue< T , float > >::const_iterator mapi = depMap->find(recoecalcands[i]);
00607             if (mapi!=depMap->end()){  // found candidate in isolation map! 
00608               etahistiso[n+1]->Fill(recoecalcands[i]->eta(),mapi->val);
00609               phihistiso[n+1]->Fill(recoecalcands[i]->phi(),mapi->val);
00610               ethistiso[n+1]->Fill(recoecalcands[i]->et(),mapi->val);
00611             }
00612           }
00613         }
00614       }
00615     } // END of if n+1 < then the number of hlt collections
00616   }
00617 
00618 
00620   //        Fill mc matched objects into histograms         //
00622   unsigned int mtachedMcParts = 0;
00623   float mindist=0.3;
00624   if(n==0) mindist=0.5; //low L1-resolution => allow wider matching 
00625   for(unsigned int i =0; i < gencut_; i++){
00626     //match generator candidate    
00627     bool matchThis= false;
00628     math::XYZVector candDir=sortedGen[i].momentum();
00629     unsigned int closest = 0;
00630     double closestDr = 1000.;
00631     for(unsigned int trigOb = 0 ; trigOb < recoecalcands.size(); trigOb++){
00632       double dr = DeltaR(recoecalcands[trigOb]->momentum(),candDir);
00633       if (dr < closestDr) {
00634         closestDr = dr;
00635         closest = trigOb;
00636       }
00637       if (closestDr > mindist) { // it's not really a "match" if it's that far away
00638         closest = -1;
00639       } else {
00640         mtachedMcParts++;
00641         matchThis = true;
00642       }
00643     }
00644     if ( !matchThis ) continue; // only plot matched candidates
00645     // fill coordinates of mc particle matching trigger object
00646     ethistmatch[n] ->Fill( sortedGen[i].et()  );
00647     if (sortedGen[i].et() > plotMinEtForEtaEffPlot) {
00648       etahistmatch[n]->Fill( sortedGen[i].eta() );
00649       phihistmatch[n]->Fill( sortedGen[i].phi() );
00650     }
00652     //  Plot isolation variables (show the not-yet-cut        //
00653     //  isolation, i.e. associated to next filter)            //
00655     if (n+1 < numOfHLTCollectionLabels){ // can't plot beyond last
00656       if (plotiso[n+1] ){  // only plot if requested in config
00657         for (unsigned int j =  0 ; j < isoNames[n+1].size() ;j++  ){
00658           edm::Handle<edm::AssociationMap<edm::OneToValue< T , float > > > depMap; 
00659           iEvent.getByLabel(isoNames[n+1].at(j),depMap);
00660           if (depMap.isValid()){ //Map may not exist if only one candidate passes a double filter
00661             typename edm::AssociationMap<edm::OneToValue< T , float > >::const_iterator mapi = depMap->find(recoecalcands[closest]);
00662             if (mapi!=depMap->end()){  // found candidate in isolation map!
00663               // Only make efficiency plot using photons with some min Et
00664               etahistisomatch[n+1]->Fill(sortedGen[i].eta(),mapi->val);
00665               phihistisomatch[n+1]->Fill(sortedGen[i].phi(),mapi->val);
00666               ethistisomatch[n+1]->Fill(sortedGen[i].et(),mapi->val);
00667             }
00668           }
00669         }
00670       }
00671     } // END of if n+1 < then the number of hlt collections
00672   }
00673   // fill total mc matched efficiency
00674   if (mtachedMcParts >= reqNum ) 
00675     totalmatch->Fill(n+0.5);
00676   
00677 
00678 }
00679 
00680 
00682 //      method called once each job just after ending the event loop          //
00684 void EmDQM::endJob(){
00685 
00686 }
00687 
00688 DEFINE_FWK_MODULE(EmDQM);