Reference Manual

Go to the documentation of this file.
00001 
00008 // Original Author:  Zoltan Gecse
00009 
00010 #include <memory>
00011 #include "FWCore/Framework/interface/Frameworkfwd.h"
00012 #include "FWCore/Framework/interface/EDAnalyzer.h"
00013 #include "FWCore/Framework/interface/Event.h"
00014 #include "FWCore/Framework/interface/MakerMacros.h"
00015 #include "FWCore/MessageLogger/interface/MessageLogger.h"
00016 #include "FWCore/ServiceRegistry/interface/Service.h"
00017 #include "FWCore/ParameterSet/interface/ParameterSet.h"
00018 #include "FWCore/Utilities/interface/InputTag.h"
00019 #include "FWCore/Common/interface/TriggerNames.h"
00020 
00021 #include "DataFormats/Common/interface/RefToBase.h"
00022 #include "DataFormats/TrackReco/interface/Track.h"
00023 #include "DataFormats/RecoCandidate/interface/RecoChargedCandidate.h"
00024 #include "DataFormats/RecoCandidate/interface/RecoChargedCandidateFwd.h"
00025 #include "DataFormats/L1Trigger/interface/L1MuonParticle.h"
00026 #include "DataFormats/L1Trigger/interface/L1MuonParticleFwd.h"
00027 #include "DataFormats/Candidate/interface/Candidate.h"
00028 #include "DataFormats/HepMCCandidate/interface/GenParticle.h"
00029 #include "DataFormats/MuonReco/interface/MuonFwd.h"
00030 #include "DataFormats/MuonReco/interface/Muon.h"
00031 #include "DataFormats/Math/interface/deltaR.h"
00032 #include "DataFormats/Common/interface/TriggerResults.h"
00033 #include "DataFormats/Candidate/interface/Particle.h"
00034 #include "DataFormats/HLTReco/interface/TriggerEventWithRefs.h"
00035 #include "DataFormats/HLTReco/interface/TriggerEvent.h"
00036 
00037 #include "HLTrigger/HLTcore/interface/HLTConfigProvider.h"
00038 
00039 #include "DQMServices/Core/interface/DQMStore.h"
00040 #include "DQMServices/Core/interface/MonitorElement.h"
00041 
00042 #include "CommonTools/Utils/interface/PtComparator.h"
00043 
00044 #include "TLorentzVector.h"
00045 
00046 using namespace std;
00047 using namespace edm;
00048 using namespace reco;
00049 using namespace l1extra;
00050 using namespace trigger;
00051 
00052 class HeavyFlavorValidation : public edm::EDAnalyzer {
00053   public:
00054     explicit HeavyFlavorValidation(const edm::ParameterSet&);
00055     ~HeavyFlavorValidation();
00056   private:
00057     virtual void beginRun(const Run & iRun, const EventSetup & iSetup);
00058     virtual void analyze(const edm::Event&, const edm::EventSetup&);
00059     virtual void endJob();
00060     int getMotherId( const Candidate * p );
00061     void match( MonitorElement * me, vector<LeafCandidate> & from, vector<LeafCandidate> & to, double deltaRMatchingCut, vector<int> & map );
00062     void myBook2D( TString name, vector<double> &xBins, TString xLabel, vector<double> &yBins, TString yLabel, TString title);
00063     void myBook2D( TString name, vector<double> &xBins, TString xLabel, vector<double> &yBins, TString yLabel){
00064       myBook2D( name, xBins, xLabel, yBins, yLabel, name);
00065     }
00066     void myBookProfile2D( TString name, vector<double> &xBins, TString xLabel, vector<double> &yBins, TString yLabel, TString title);
00067     void myBookProfile2D( TString name, vector<double> &xBins, TString xLabel, vector<double> &yBins, TString yLabel){
00068       myBookProfile2D( name, xBins, xLabel, yBins, yLabel, name);
00069     }
00070     void myBook1D( TString name, vector<double> &xBins, TString label, TString title );
00071     void myBook1D( TString name, vector<double> &xBins, TString label ){
00072       myBook1D( name, xBins, label, name );
00073     }
00074     string dqmFolder;
00075     string triggerProcessName;
00076     string triggerPathName;
00077     InputTag triggerSummaryRAWTag;
00078     InputTag triggerSummaryAODTag;
00079     InputTag triggerResultsTag;
00080     InputTag recoMuonsTag;
00081     InputTag genParticlesTag;   
00082     vector<int> motherIDs;
00083     double genGlobDeltaRMatchingCut;
00084     double globL1DeltaRMatchingCut;
00085     double globL2DeltaRMatchingCut;
00086     double globL3DeltaRMatchingCut;
00087     vector<double> deltaEtaBins;
00088     vector<double> deltaPhiBins;
00089     vector<double> muonPtBins;
00090     vector<double> muonEtaBins;
00091     vector<double> muonPhiBins;
00092     vector<double> dimuonPtBins;
00093     vector<double> dimuonEtaBins;
00094     vector<double> dimuonDRBins;
00095     DQMStore* dqmStore;
00096     map<TString, MonitorElement *> ME;
00097     vector<pair<string,int> > filterNamesLevels;
00098     const double muonMass;
00099 };
00100 
00101 HeavyFlavorValidation::HeavyFlavorValidation(const ParameterSet& pset):
00102 //get parameters
00103   dqmFolder(pset.getUntrackedParameter<string>("DQMFolder")),
00104   triggerProcessName(pset.getUntrackedParameter<string>("TriggerProcessName")),
00105   triggerPathName(pset.getUntrackedParameter<string>("TriggerPathName")),
00106   triggerSummaryRAWTag(InputTag( pset.getUntrackedParameter<string>("TriggerSummaryRAW"), "", triggerProcessName)),
00107   triggerSummaryAODTag(InputTag( pset.getUntrackedParameter<string>("TriggerSummaryAOD"), "", triggerProcessName)),
00108   triggerResultsTag(InputTag( pset.getUntrackedParameter<string>("TriggerResults"), "", triggerProcessName)),
00109   recoMuonsTag(pset.getParameter<InputTag>("RecoMuons")),
00110   genParticlesTag(pset.getParameter<InputTag>("GenParticles")),
00111   motherIDs(pset.getUntrackedParameter<vector<int> >("MotherIDs")),
00112   genGlobDeltaRMatchingCut(pset.getUntrackedParameter<double>("GenGlobDeltaRMatchingCut")),
00113   globL1DeltaRMatchingCut(pset.getUntrackedParameter<double>("GlobL1DeltaRMatchingCut")),
00114   globL2DeltaRMatchingCut(pset.getUntrackedParameter<double>("GlobL2DeltaRMatchingCut")),
00115   globL3DeltaRMatchingCut(pset.getUntrackedParameter<double>("GlobL3DeltaRMatchingCut")),
00116   deltaEtaBins(pset.getUntrackedParameter<vector<double> >("DeltaEtaBins")),
00117   deltaPhiBins(pset.getUntrackedParameter<vector<double> >("DeltaPhiBins")),
00118   muonPtBins(pset.getUntrackedParameter<vector<double> >("MuonPtBins")),
00119   muonEtaBins(pset.getUntrackedParameter<vector<double> >("MuonEtaBins")),
00120   muonPhiBins(pset.getUntrackedParameter<vector<double> >("MuonPhiBins")),
00121   dimuonPtBins(pset.getUntrackedParameter<vector<double> >("DimuonPtBins")),
00122   dimuonEtaBins(pset.getUntrackedParameter<vector<double> >("DimuonEtaBins")),
00123   dimuonDRBins(pset.getUntrackedParameter<vector<double> >("DimuonDRBins")),
00124   muonMass(0.106)
00125 {}
00126   
00127 void HeavyFlavorValidation::beginRun(const Run & iRun, const EventSetup & iSetup){
00128 //discover HLT configuration
00129   HLTConfigProvider hltConfig;
00130   bool isChanged;
00131   if(hltConfig.init(iRun, iSetup, triggerProcessName, isChanged)){
00132     LogDebug("HLTriggerOfflineHeavyFlavor") << "Successfully initialized HLTConfigProvider with process name: "<<triggerProcessName<<endl;
00133   }else{
00134     LogWarning("HLTriggerOfflineHeavyFlavor") << "Could not initialize HLTConfigProvider with process name: "<<triggerProcessName<<endl;
00135     return;
00136   }
00137   stringstream os;
00138   vector<string> triggerNames = hltConfig.triggerNames();
00139   for( size_t i = 0; i < triggerNames.size(); i++) {
00140     TString triggerName = triggerNames[i];
00141     if (triggerName.Contains(triggerPathName)){ 
00142       vector<string> moduleNames = hltConfig.moduleLabels( triggerNames[i] );
00143       for( size_t j = 0; j < moduleNames.size(); j++) {
00144         TString name = moduleNames[j];
00145         if(name.Contains("Filtered")){
00146           int level = 0;
00147           if(name.Contains("L1"))
00148             level = 1;
00149           else if(name.Contains("L2"))
00150             level = 2;
00151           else if(name.Contains("L3"))
00152             level = 3;
00153           filterNamesLevels.push_back( pair<string,int>(moduleNames[j],level) );
00154           os<<" "<<moduleNames[j];
00155         }
00156       }
00157       break;
00158     }
00159   }
00160   if(filterNamesLevels.size()==0){
00161     LogDebug("HLTriggerOfflineHeavyFlavor")<<"Bad Trigger Path: "<<triggerPathName<<endl;
00162     return;
00163   }else{
00164     LogDebug("HLTriggerOfflineHeavyFlavor")<<"Trigger Path: "<<triggerPathName<<" has filters:"<<os.str();
00165   }
00166 
00167 //create Monitor Elements
00168   dqmStore = Service<DQMStore>().operator->();  
00169   if( !dqmStore ){
00170     LogError("HLTriggerOfflineHeavyFlavor") << "Could not find DQMStore service\n";
00171     return;
00172   }
00173   dqmStore->setVerbose(0);
00174   dqmStore->setCurrentFolder((dqmFolder+"/")+triggerProcessName+"/"+triggerPathName);
00175 // Eta Pt Single  
00176   myBook2D( "genMuon_genEtaPt", muonEtaBins, "#mu eta", muonPtBins, " #mu pT (GeV)");
00177   myBook2D( "globMuon_genEtaPt", muonEtaBins, "#mu eta", muonPtBins, " #mu pT (GeV)");
00178   myBook2D( "globMuon_recoEtaPt", muonEtaBins, "#mu eta", muonPtBins, " #mu pT (GeV)");
00179   for(size_t i=0; i<filterNamesLevels.size(); i++){
00180     myBook2D( TString::Format("filt%dMuon_recoEtaPt",int(i+1)), muonEtaBins, "#mu eta", muonPtBins, " #mu pT (GeV)", filterNamesLevels[i].first);
00181   }
00182   myBook2D( "pathMuon_recoEtaPt", muonEtaBins, "#mu eta", muonPtBins, " #mu pT (GeV)", triggerPathName);
00183   myBook2D( "resultMuon_recoEtaPt", muonEtaBins, "#mu eta", muonPtBins, " #mu pT (GeV)");
00184 // Eta Pt Single Resolution
00185   myBookProfile2D( "resGlobGen_genEtaPt", muonEtaBins, "#mu eta", muonPtBins, " #mu pT (GeV)");
00186   for(size_t i=0; i<filterNamesLevels.size(); i++){
00187     myBookProfile2D( TString::Format("resFilt%dGlob_recoEtaPt",int(i+1)), muonEtaBins, "#mu eta", muonPtBins, " #mu pT (GeV)", filterNamesLevels[i].first);
00188   }
00189   myBookProfile2D( "resPathGlob_recoEtaPt", muonEtaBins, "#mu eta", muonPtBins, " #mu pT (GeV)", triggerPathName);
00190 // Eta Pt Double  
00191   myBook2D( "genDimuon_genEtaPt", dimuonEtaBins, "#mu#mu eta", dimuonPtBins, " #mu#mu pT (GeV)");
00192   myBook2D( "globDimuon_genEtaPt", dimuonEtaBins, "#mu#mu eta", dimuonPtBins, " #mu#mu pT (GeV)");
00193   myBook2D( "globDimuon_recoEtaPt", dimuonEtaBins, "#mu#mu eta", dimuonPtBins, " #mu#mu pT (GeV)");
00194   for(size_t i=0; i<filterNamesLevels.size(); i++){
00195     myBook2D( TString::Format("filt%dDimuon_recoEtaPt",int(i+1)), dimuonEtaBins, "#mu#mu eta", dimuonPtBins, " #mu#mu pT (GeV)", filterNamesLevels[i].first);
00196   }  
00197   myBook2D( "pathDimuon_recoEtaPt", dimuonEtaBins, "#mu#mu eta", dimuonPtBins, " #mu#mu pT (GeV)", triggerPathName);
00198   myBook2D( "resultDimuon_recoEtaPt", dimuonEtaBins, "#mu#mu eta", dimuonPtBins, " #mu#mu pT (GeV)");
00199   for(size_t i=0; i<filterNamesLevels.size(); i++){
00200     myBook2D( TString::Format("diFilt%dDimuon_recoEtaPt",int(i+1)), dimuonEtaBins, "#mu#mu eta", dimuonPtBins, " #mu#mu pT (GeV)", filterNamesLevels[i].first);
00201   }  
00202   myBook2D( "diPathDimuon_recoEtaPt", dimuonEtaBins, "#mu#mu eta", dimuonPtBins, " #mu#mu pT (GeV)", triggerPathName);
00203 // Eta Phi Single
00204   myBook2D( "genMuon_genEtaPhi", muonEtaBins, "#mu eta", muonPhiBins, "#mu phi");
00205   myBook2D( "globMuon_genEtaPhi", muonEtaBins, "#mu eta", muonPhiBins, "#mu phi");
00206   myBook2D( "globMuon_recoEtaPhi", muonEtaBins, "#mu eta", muonPhiBins, "#mu phi");
00207   for(size_t i=0; i<filterNamesLevels.size(); i++){
00208     myBook2D( TString::Format("filt%dMuon_recoEtaPhi",int(i+1)), muonEtaBins, "#mu eta", muonPhiBins, "#mu phi", filterNamesLevels[i].first);
00209   }
00210   myBook2D( "pathMuon_recoEtaPhi", muonEtaBins, "#mu eta", muonPhiBins, "#mu phi", triggerPathName);
00211   myBook2D( "resultMuon_recoEtaPhi", muonEtaBins, "#mu eta", muonPhiBins, "#mu phi");
00212 // Rap Pt Double
00213   myBook2D( "genDimuon_genRapPt", dimuonEtaBins, "#mu#mu rapidity", dimuonPtBins, " #mu#mu pT (GeV)");
00214   myBook2D( "globDimuon_genRapPt", dimuonEtaBins, "#mu#mu rapidity", dimuonPtBins, " #mu#mu pT (GeV)");
00215   myBook2D( "globDimuon_recoRapPt", dimuonEtaBins, "#mu#mu rapidity", dimuonPtBins, " #mu#mu pT (GeV)");
00216   for(size_t i=0; i<filterNamesLevels.size(); i++){
00217     myBook2D( TString::Format("filt%dDimuon_recoRapPt",int(i+1)), dimuonEtaBins, "#mu#mu rapidity", dimuonPtBins, " #mu#mu pT (GeV)", filterNamesLevels[i].first);
00218   }  
00219   myBook2D( "pathDimuon_recoRapPt", dimuonEtaBins, "#mu#mu rapidity", dimuonPtBins, " #mu#mu pT (GeV)", triggerPathName);
00220   myBook2D( "resultDimuon_recoRapPt", dimuonEtaBins, "#mu#mu rapidity", dimuonPtBins, " #mu#mu pT (GeV)");
00221   for(size_t i=0; i<filterNamesLevels.size(); i++){
00222     myBook2D( TString::Format("diFilt%dDimuon_recoRapPt",int(i+1)), dimuonEtaBins, "#mu#mu rapidity", dimuonPtBins, " #mu#mu pT (GeV)", filterNamesLevels[i].first);
00223   }  
00224   myBook2D( "diPathDimuon_recoRapPt", dimuonEtaBins, "#mu#mu rapidity", dimuonPtBins, " #mu#mu pT (GeV)", triggerPathName);
00225 // Pt DR Double  
00226   myBook2D( "genDimuon_genPtDR", dimuonPtBins, " #mu#mu pT (GeV)", dimuonDRBins, "#mu#mu #Delta R at IP");
00227   myBook2D( "globDimuon_genPtDR", dimuonPtBins, " #mu#mu pT (GeV)", dimuonDRBins, "#mu#mu #Delta R at IP");
00228   myBook2D( "globDimuon_recoPtDR", dimuonPtBins, " #mu#mu pT (GeV)", dimuonDRBins, "#mu#mu #Delta R at IP");
00229   for(size_t i=0; i<filterNamesLevels.size(); i++){
00230     myBook2D( TString::Format("filt%dDimuon_recoPtDR",int(i+1)), dimuonPtBins, " #mu#mu pT (GeV)", dimuonDRBins, "#mu#mu #Delta R at IP", filterNamesLevels[i].first);
00231   }
00232   myBook2D( "pathDimuon_recoPtDR", dimuonPtBins, " #mu#mu pT (GeV)", dimuonDRBins, "#mu#mu #Delta R at IP", triggerPathName);
00233   for(size_t i=0; i<filterNamesLevels.size(); i++){
00234     myBook2D( TString::Format("diFilt%dDimuon_recoPtDR",int(i+1)), dimuonPtBins, " #mu#mu pT (GeV)", dimuonDRBins, "#mu#mu #Delta R at IP", filterNamesLevels[i].first);
00235   }  
00236   myBook2D( "diPathDimuon_recoPtDR", dimuonPtBins, " #mu#mu pT (GeV)", dimuonDRBins, "#mu#mu #Delta R at IP", triggerPathName);
00237   myBook2D( "resultDimuon_recoPtDR", dimuonPtBins, " #mu#mu pT (GeV)", dimuonDRBins, "#mu#mu #Delta R at IP");
00238 // Pt DRpos Double  
00239 //   myBook2D( "genDimuon_genPtDRpos", dimuonPtBins, " #mu#mu pT (GeV)", dimuonDRBins, "#mu#mu #Delta R in MS");
00240 //   myBook2D( "globDimuon_genPtDRpos", dimuonPtBins, " #mu#mu pT (GeV)", dimuonDRBins, "#mu#mu #Delta R in MS");
00241   myBook2D( "globDimuon_recoPtDRpos", dimuonPtBins, " #mu#mu pT (GeV)", dimuonDRBins, "#mu#mu #Delta R in MS");
00242   for(size_t i=0; i<filterNamesLevels.size(); i++){
00243     myBook2D( TString::Format("filt%dDimuon_recoPtDRpos",int(i+1)), dimuonPtBins, " #mu#mu pT (GeV)", dimuonDRBins, "#mu#mu #Delta R in MS", filterNamesLevels[i].first);
00244   }
00245   myBook2D( "pathDimuon_recoPtDRpos", dimuonPtBins, " #mu#mu pT (GeV)", dimuonDRBins, "#mu#mu #Delta R in MS", triggerPathName);
00246   for(size_t i=0; i<filterNamesLevels.size(); i++){
00247     myBook2D( TString::Format("diFilt%dDimuon_recoPtDRpos",int(i+1)), dimuonPtBins, " #mu#mu pT (GeV)", dimuonDRBins, "#mu#mu #Delta R in MS", filterNamesLevels[i].first);
00248   }  
00249   myBook2D( "diPathDimuon_recoPtDRpos", dimuonPtBins, " #mu#mu pT (GeV)", dimuonDRBins, "#mu#mu #Delta R in MS", triggerPathName);
00250   myBook2D( "resultDimuon_recoPtDRpos", dimuonPtBins, " #mu#mu pT (GeV)", dimuonDRBins, "#mu#mu #Delta R in MS");
00251 
00252 /*  ME["genDimuon_genPt+Pt-"] = myBook2D( "genDimuon_genPt+Pt-", muonPtBins, "#mu+ pt (GeV)", muonPtBins, "#mu- pT (GeV)");
00253   ME["genGlobDimuon_genPt+Pt-"] = myBook2D( "genGlobDimuon_genPt+Pt-", muonPtBins, "#mu+ pt (GeV)", muonPtBins, "#mu- pT (GeV)");
00254   ME["genGlobDimuon_recoPt+Pt-"] = myBook2D( "genGlobDimuon_recoPt+Pt-", muonPtBins, "#mu+ pt (GeV)", muonPtBins, "#mu- pT (GeV)");
00255   ME["genGlobL1Dimuon_recoPt+Pt-"] = myBook2D( "genGlobL1Dimuon_recoPt+Pt-", muonPtBins, "#mu+ pt (GeV)", muonPtBins, "#mu- pT (GeV)");
00256   ME["genGlobL1L2Dimuon_recoPt+Pt-"] = myBook2D( "genGlobL1L2Dimuon_recoPt+Pt-", muonPtBins, "#mu+ pt (GeV)", muonPtBins, "#mu- pT (GeV)");
00257   ME["genGlobL1L2L3Dimuon_recoPt+Pt-"] = myBook2D( "genGlobL1L2L3Dimuon_recoPt+Pt-", muonPtBins, "#mu+ pt (GeV)", muonPtBins, "#mu- pT (GeV)");*/
00258 
00259 // Matching
00260   myBook2D( "globGen_deltaEtaDeltaPhi", deltaEtaBins, "#Delta eta", deltaPhiBins, "#Delta phi");
00261   for(size_t i=0; i<filterNamesLevels.size(); i++){
00262     myBook2D( TString::Format("filt%dGlob_deltaEtaDeltaPhi",int(i+1)), deltaEtaBins, "#Delta eta", deltaPhiBins, "#Delta phi", filterNamesLevels[i].first);
00263   }    
00264   myBook2D( "pathGlob_deltaEtaDeltaPhi", deltaEtaBins, "#Delta eta", deltaPhiBins, "#Delta phi", triggerPathName);
00265 // Size of containers
00266   vector<double> sizeBins; sizeBins.push_back(10); sizeBins.push_back(0); sizeBins.push_back(10);
00267   myBook1D( "genMuon_size", sizeBins, "container size" );
00268   myBook1D( "globMuon_size", sizeBins, "container size" );
00269   for(size_t i=0; i<filterNamesLevels.size(); i++){
00270     myBook1D( TString::Format("filt%dMuon_size",int(i+1)), sizeBins, "container size", filterNamesLevels[i].first);
00271   }    
00272   myBook1D( "pathMuon_size", sizeBins, "container size", triggerPathName );
00273 }
00274 
00275 void HeavyFlavorValidation::analyze(const Event& iEvent, const EventSetup& iSetup){
00276   if( !dqmStore ){
00277     LogDebug("HLTriggerOfflineHeavyFlavor")<<"Could not access DQM Store service"<<endl;
00278     return;
00279   }
00280   if(filterNamesLevels.size()==0){
00281     return;
00282   }
00283 //access the containers and create LeafCandidate copies
00284   vector<LeafCandidate> genMuons;
00285   Handle<GenParticleCollection> genParticles;
00286   iEvent.getByLabel(genParticlesTag, genParticles);
00287   if(genParticles.isValid()){
00288     for(GenParticleCollection::const_iterator p=genParticles->begin(); p!= genParticles->end(); ++p){
00289       if( p->status() == 1 && std::abs(p->pdgId())==13 && 
00290           ( find( motherIDs.begin(), motherIDs.end(), -1 )!=motherIDs.end() || find( motherIDs.begin(), motherIDs.end(), getMotherId(&(*p)) )!=motherIDs.end() ) ){
00291         genMuons.push_back( *p );
00292       }  
00293     }
00294   }else{
00295     LogDebug("HLTriggerOfflineHeavyFlavor")<<"Could not access GenParticleCollection"<<endl;
00296   }
00297   sort(genMuons.begin(), genMuons.end(), GreaterByPt<LeafCandidate>());
00298   ME["genMuon_size"]->Fill(genMuons.size());
00299   LogDebug("HLTriggerOfflineHeavyFlavor")<<"GenParticleCollection from "<<genParticlesTag<<" has size: "<<genMuons.size()<<endl;
00300   
00301   vector<LeafCandidate> globMuons;
00302   vector<LeafCandidate> globMuons_position;
00303   Handle<MuonCollection> recoMuonsHandle;
00304   iEvent.getByLabel(recoMuonsTag, recoMuonsHandle);
00305   if(recoMuonsHandle.isValid()){
00306     for(MuonCollection::const_iterator p=recoMuonsHandle->begin(); p!= recoMuonsHandle->end(); ++p){
00307       if(p->isGlobalMuon()){
00308         globMuons.push_back(*p);
00309         globMuons_position.push_back( LeafCandidate( p->charge(), math::XYZTLorentzVector(p->outerTrack()->innerPosition().x(), p->outerTrack()->innerPosition().y(), p->outerTrack()->innerPosition().z(), 0.) ) );
00310       }
00311     }
00312   }else{
00313     LogDebug("HLTriggerOfflineHeavyFlavor")<<"Could not access reco Muons"<<endl;
00314   }
00315   ME["globMuon_size"]->Fill(globMuons.size());
00316   LogDebug("HLTriggerOfflineHeavyFlavor")<<"Global Muons from "<<recoMuonsTag<<" has size: "<<globMuons.size()<<endl;
00317 
00318 // access RAW trigger event  
00319   vector<vector<LeafCandidate> > muonsAtFilter;
00320   vector<vector<LeafCandidate> > muonPositionsAtFilter;  
00321   for(size_t i=0; i<filterNamesLevels.size(); i++){
00322     muonsAtFilter.push_back(vector<LeafCandidate>());
00323     muonPositionsAtFilter.push_back(vector<LeafCandidate>());
00324   }
00325   Handle<TriggerEventWithRefs> rawTriggerEvent;
00326   iEvent.getByLabel( triggerSummaryRAWTag, rawTriggerEvent );
00327   if( rawTriggerEvent.isValid() ){
00328     for(size_t i=0; i<filterNamesLevels.size(); i++){
00329       size_t index = rawTriggerEvent->filterIndex(InputTag( filterNamesLevels[i].first, "", triggerProcessName ));
00330       if ( index < rawTriggerEvent->size() ){
00331         if( filterNamesLevels[i].second==1 ){
00332           vector<L1MuonParticleRef> l1Cands;
00333           rawTriggerEvent->getObjects( index, TriggerL1Mu, l1Cands );
00334           for(size_t j=0; j<l1Cands.size(); j++){
00335             muonsAtFilter[i].push_back(*l1Cands[j]);
00336           }
00337         }else{
00338           vector<RecoChargedCandidateRef> hltCands;        
00339           rawTriggerEvent->getObjects( index, TriggerMuon, hltCands );
00340           for(size_t j=0; j<hltCands.size(); j++){
00341             muonsAtFilter[i].push_back(*hltCands[j]);
00342             if( filterNamesLevels[i].second==2 ){
00343               muonPositionsAtFilter[i].push_back( LeafCandidate( hltCands[j]->charge(), math::XYZTLorentzVector(hltCands[j]->track()->innerPosition().x(), hltCands[j]->track()->innerPosition().y(), hltCands[j]->track()->innerPosition().z(), 0.) ) );
00344             }
00345           }
00346         }
00347       }
00348       ME[TString::Format("filt%dMuon_size",int(i+1))]->Fill(muonsAtFilter[i].size());
00349       LogDebug("HLTriggerOfflineHeavyFlavor")<<"Filter \""<<filterNamesLevels[i].first<<"\" has "<<muonsAtFilter[i].size()<<" muons"<<endl;
00350     }
00351   }else{
00352     LogDebug("HLTriggerOfflineHeavyFlavor")<<"Could not access RAWTriggerEvent"<<endl;
00353   }
00354 
00355 // access AOD trigger event  
00356   vector<LeafCandidate> pathMuons;
00357   Handle<TriggerEvent> aodTriggerEvent;
00358   iEvent.getByLabel(triggerSummaryAODTag,aodTriggerEvent);
00359   if(aodTriggerEvent.isValid()){
00360     TriggerObjectCollection allObjects = aodTriggerEvent->getObjects();
00361     for(int i=0; i<aodTriggerEvent->sizeFilters(); i++){         
00362      if(aodTriggerEvent->filterTag(i)==InputTag((filterNamesLevels.end()-1)->first,"",triggerProcessName)){
00363         Keys keys = aodTriggerEvent->filterKeys(i);
00364         for(size_t j=0; j<keys.size(); j++){
00365           pathMuons.push_back( LeafCandidate( allObjects[keys[j]].id()>0 ? 1 : -1, math::PtEtaPhiMLorentzVector( allObjects[keys[j]].pt(), allObjects[keys[j]].eta(), allObjects[keys[j]].phi(), muonMass ) ) );
00366         }
00367       }
00368     }
00369     ME["pathMuon_size"]->Fill(pathMuons.size());
00370     LogDebug("HLTriggerOfflineHeavyFlavor")<<"Path \""<<triggerPathName<<"\" has "<<pathMuons.size()<<" muons at last filter \""<<(filterNamesLevels.end()-1)->first<<"\""<<endl;
00371   }else{
00372     LogDebug("HLTriggerOfflineHeavyFlavor")<<"Could not access AODTriggerEvent"<<endl;
00373   }
00374 
00375 // access Trigger Results
00376   bool triggerFired = false;
00377   Handle<TriggerResults> triggerResults;
00378   iEvent.getByLabel(triggerResultsTag,triggerResults);
00379   if(triggerResults.isValid()){
00380     LogDebug("HLTriggerOfflineHeavyFlavor")<<"Successfully initialized "<<triggerResultsTag<<endl;
00381     const edm::TriggerNames & triggerNames = iEvent.triggerNames(*triggerResults);
00382     size_t index = triggerNames.triggerIndex(triggerPathName);
00383     if( index < triggerNames.size() ){
00384       triggerFired = triggerResults->accept( index );
00385     }else{
00386       LogDebug("HLTriggerOfflineHeavyFlavor")<<triggerResultsTag<<" has no trigger: "<<triggerPathName<<endl;
00387     }
00388   }else{
00389     LogDebug("HLTriggerOfflineHeavyFlavor")<<"Could not initialize "<<triggerResultsTag<<endl;
00390   }
00391 
00392 //create matching maps
00393   vector<int> glob_gen(genMuons.size(),-1);
00394   match( ME["globGen_deltaEtaDeltaPhi"], genMuons, globMuons, genGlobDeltaRMatchingCut, glob_gen );
00395   vector<vector<int> > filt_glob;
00396   for(size_t i=0; i<filterNamesLevels.size(); i++){
00397     filt_glob.push_back( vector<int>(globMuons.size(),-1) );            
00398     if( filterNamesLevels[i].second == 1 ){
00399       match( ME[TString::Format("filt%dGlob_deltaEtaDeltaPhi",int(i+1))], globMuons_position, muonsAtFilter[i] ,globL1DeltaRMatchingCut, filt_glob[i] );
00400     }else if( filterNamesLevels[i].second == 2 ){
00401       match( ME[TString::Format("filt%dGlob_deltaEtaDeltaPhi",int(i+1))], globMuons_position, muonPositionsAtFilter[i] ,globL2DeltaRMatchingCut, filt_glob[i] );
00402     }else if( filterNamesLevels[i].second == 3 ){
00403       match( ME[TString::Format("filt%dGlob_deltaEtaDeltaPhi",int(i+1))], globMuons, muonsAtFilter[i] ,globL3DeltaRMatchingCut, filt_glob[i] );
00404     }
00405   }
00406   vector<int> path_glob(globMuons.size(),-1);
00407   if( (filterNamesLevels.end()-1)->second == 1 ){
00408     match( ME["pathGlob_deltaEtaDeltaPhi"], globMuons_position, pathMuons ,globL1DeltaRMatchingCut, path_glob );
00409   }else if( (filterNamesLevels.end()-1)->second == 2 ){
00410     match( ME["pathGlob_deltaEtaDeltaPhi"], globMuons, pathMuons ,globL2DeltaRMatchingCut, path_glob );
00411   }else if( (filterNamesLevels.end()-1)->second == 3 ){
00412     match( ME["pathGlob_deltaEtaDeltaPhi"], globMuons, pathMuons ,globL3DeltaRMatchingCut, path_glob );
00413   }
00414     
00415 //fill histos
00416   bool first = true;
00417   for(size_t i=0; i<genMuons.size(); i++){
00418     ME["genMuon_genEtaPt"]->Fill(genMuons[i].eta(), genMuons[i].pt());
00419     ME["genMuon_genEtaPhi"]->Fill(genMuons[i].eta(), genMuons[i].phi());
00420     if(glob_gen[i] != -1){
00421       ME["resGlobGen_genEtaPt"]->Fill(genMuons[i].eta(), genMuons[i].pt(), (globMuons[glob_gen[i]].pt()-genMuons[i].pt())/genMuons[i].pt() );
00422       ME["globMuon_genEtaPt"]->Fill(genMuons[i].eta(), genMuons[i].pt());
00423       ME["globMuon_genEtaPhi"]->Fill(genMuons[i].eta(), genMuons[i].phi());
00424       ME["globMuon_recoEtaPt"]->Fill(globMuons[glob_gen[i]].eta(), globMuons[glob_gen[i]].pt());
00425       ME["globMuon_recoEtaPhi"]->Fill(globMuons[glob_gen[i]].eta(), globMuons[glob_gen[i]].phi());
00426       for(size_t f=0; f<filterNamesLevels.size(); f++){
00427         if(filt_glob[f][glob_gen[i]] != -1){
00428           ME[TString::Format("resFilt%dGlob_recoEtaPt",int(f+1))]->Fill(globMuons[glob_gen[i]].eta(), globMuons[glob_gen[i]].pt(), (muonsAtFilter[f][filt_glob[f][glob_gen[i]]].pt()-globMuons[glob_gen[i]].pt())/globMuons[glob_gen[i]].pt() );
00429           ME[TString::Format("filt%dMuon_recoEtaPt",int(f+1))]->Fill(globMuons[glob_gen[i]].eta(), globMuons[glob_gen[i]].pt());
00430           ME[TString::Format("filt%dMuon_recoEtaPhi",int(f+1))]->Fill(globMuons[glob_gen[i]].eta(), globMuons[glob_gen[i]].phi());
00431         }else{
00432           break;
00433         }  
00434       }
00435       if(path_glob[glob_gen[i]] != -1){
00436         ME["resPathGlob_recoEtaPt"]->Fill(globMuons[glob_gen[i]].eta(), globMuons[glob_gen[i]].pt(), (pathMuons[path_glob[glob_gen[i]]].pt()-globMuons[glob_gen[i]].pt())/globMuons[glob_gen[i]].pt() );
00437         ME["pathMuon_recoEtaPt"]->Fill(globMuons[glob_gen[i]].eta(), globMuons[glob_gen[i]].pt());
00438         ME["pathMuon_recoEtaPhi"]->Fill(globMuons[glob_gen[i]].eta(), globMuons[glob_gen[i]].phi());
00439       }
00440 //highest pt muon      
00441       if( first ){
00442         first = false;
00443         if( triggerFired ){
00444           ME["resultMuon_recoEtaPt"]->Fill(globMuons[glob_gen[i]].eta(), globMuons[glob_gen[i]].pt());
00445           ME["resultMuon_recoEtaPhi"]->Fill(globMuons[glob_gen[i]].eta(), globMuons[glob_gen[i]].phi());
00446         }
00447       }
00448     }
00449   }  
00450 
00451 //fill dimuon histograms (highest pT, opposite charge) 
00452   int secondMuon = 0;
00453   for(size_t j=1; j<genMuons.size(); j++){
00454     if(genMuons[0].charge()*genMuons[j].charge()==-1){
00455       secondMuon = j;
00456       break;
00457     }
00458   }
00459   if(secondMuon > 0){
00460 //    int pos = genMuons[0].charge()>0 ? 0 : secondMuon ;
00461 //    int neg = genMuons[0].charge()<0 ? 0 : secondMuon ;
00462 //two generated
00463     double genDimuonPt = (genMuons[0].p4()+genMuons[secondMuon].p4()).pt();
00464     double genDimuonEta = (genMuons[0].p4()+genMuons[secondMuon].p4()).eta();
00465     double genDimuonRap = (genMuons[0].p4()+genMuons[secondMuon].p4()).Rapidity();
00466     double genDimuonDR = deltaR<LeafCandidate,LeafCandidate>( genMuons[0], genMuons[secondMuon] );
00467     bool highPt = genMuons[0].pt()>7. && genMuons[secondMuon].pt()>7;
00468     ME["genDimuon_genEtaPt"]->Fill( genDimuonEta, genDimuonPt );
00469     ME["genDimuon_genRapPt"]->Fill( genDimuonRap, genDimuonPt );
00470     if(highPt) ME["genDimuon_genPtDR"]->Fill( genDimuonPt, genDimuonDR );
00471 //two global
00472     if(glob_gen[0]!=-1 && glob_gen[secondMuon]!=-1){
00473       ME["globDimuon_genEtaPt"]->Fill( genDimuonEta, genDimuonPt );
00474       ME["globDimuon_genRapPt"]->Fill( genDimuonRap, genDimuonPt );
00475       if(highPt) ME["globDimuon_genPtDR"]->Fill( genDimuonPt, genDimuonDR );
00476       double globDimuonPt = (globMuons[glob_gen[0]].p4()+globMuons[glob_gen[secondMuon]].p4()).pt();
00477       double globDimuonEta = (globMuons[glob_gen[0]].p4()+globMuons[glob_gen[secondMuon]].p4()).eta();
00478       double globDimuonRap = (globMuons[glob_gen[0]].p4()+globMuons[glob_gen[secondMuon]].p4()).Rapidity();
00479       double globDimuonDR = deltaR<LeafCandidate,LeafCandidate>( globMuons[glob_gen[0]], globMuons[glob_gen[secondMuon]] );
00480       double globDimuonDRpos = deltaR<LeafCandidate,LeafCandidate>( globMuons_position[glob_gen[0]], globMuons_position[glob_gen[secondMuon]] );
00481       ME["globDimuon_recoEtaPt"]->Fill( globDimuonEta, globDimuonPt );
00482       ME["globDimuon_recoRapPt"]->Fill( globDimuonRap, globDimuonPt );
00483       if(highPt) ME["globDimuon_recoPtDR"]->Fill( globDimuonPt, globDimuonDR );
00484       if(highPt) ME["globDimuon_recoPtDRpos"]->Fill( globDimuonPt, globDimuonDRpos );
00485 //two filter objects    
00486       for(size_t f=0; f<filterNamesLevels.size(); f++){
00487         if(filt_glob[f][glob_gen[0]] != -1 && filt_glob[f][glob_gen[secondMuon]] != -1){
00488           ME[TString::Format("diFilt%dDimuon_recoEtaPt",int(f+1))]->Fill( globDimuonEta, globDimuonPt );
00489           ME[TString::Format("diFilt%dDimuon_recoRapPt",int(f+1))]->Fill( globDimuonRap, globDimuonPt );
00490           if(highPt) ME[TString::Format("diFilt%dDimuon_recoPtDR",int(f+1))]->Fill( globDimuonPt, globDimuonDR );
00491           if(highPt) ME[TString::Format("diFilt%dDimuon_recoPtDRpos",int(f+1))]->Fill( globDimuonPt, globDimuonDRpos );
00492         }else{
00493           break;
00494         }  
00495       } 
00496 //one filter object
00497       for(size_t f=0; f<filterNamesLevels.size(); f++){
00498         if(filt_glob[f][glob_gen[0]] != -1 || filt_glob[f][glob_gen[secondMuon]] != -1){
00499           ME[TString::Format("filt%dDimuon_recoEtaPt",int(f+1))]->Fill( globDimuonEta, globDimuonPt );
00500           ME[TString::Format("filt%dDimuon_recoRapPt",int(f+1))]->Fill( globDimuonRap, globDimuonPt );
00501           if(highPt) ME[TString::Format("filt%dDimuon_recoPtDR",int(f+1))]->Fill( globDimuonPt, globDimuonDR );
00502           if(highPt) ME[TString::Format("filt%dDimuon_recoPtDRpos",int(f+1))]->Fill( globDimuonPt, globDimuonDRpos );
00503         }else{
00504           break;
00505         }
00506       }
00507 //two path objects
00508       if(path_glob[glob_gen[0]] != -1 && path_glob[glob_gen[secondMuon]] != -1){
00509         ME["diPathDimuon_recoEtaPt"]->Fill( globDimuonEta, globDimuonPt );
00510         ME["diPathDimuon_recoRapPt"]->Fill( globDimuonRap, globDimuonPt );
00511         if(highPt) ME["diPathDimuon_recoPtDR"]->Fill( globDimuonPt, globDimuonDR );
00512         if(highPt) ME["diPathDimuon_recoPtDRpos"]->Fill( globDimuonPt, globDimuonDRpos );
00513       }
00514 //one path object
00515       if(path_glob[glob_gen[0]] != -1 || path_glob[glob_gen[secondMuon]] != -1){
00516         ME["pathDimuon_recoEtaPt"]->Fill( globDimuonEta, globDimuonPt );
00517         ME["pathDimuon_recoRapPt"]->Fill( globDimuonRap, globDimuonPt );
00518         if(highPt) ME["pathDimuon_recoPtDR"]->Fill( globDimuonPt, globDimuonDR );
00519         if(highPt) ME["pathDimuon_recoPtDRpos"]->Fill( globDimuonPt, globDimuonDRpos );
00520       }
00521 //trigger result
00522       if( triggerFired ){
00523         ME["resultDimuon_recoEtaPt"]->Fill( globDimuonEta, globDimuonPt );
00524         ME["resultDimuon_recoRapPt"]->Fill( globDimuonRap, globDimuonPt );
00525         if(highPt) ME["resultDimuon_recoPtDR"]->Fill( globDimuonPt, globDimuonDR );
00526         if(highPt) ME["resultDimuon_recoPtDRpos"]->Fill( globDimuonPt, globDimuonDRpos );
00527       }
00528     }
00529   }
00530 }
00531 
00532 void HeavyFlavorValidation::endJob(){
00533 }
00534 
00535 int HeavyFlavorValidation::getMotherId( const Candidate * p ){
00536   const Candidate* mother = p->mother();
00537   if( mother ){
00538     if( mother->pdgId() == p->pdgId() ){
00539       return getMotherId(mother);
00540     }else{
00541       return mother->pdgId();
00542     }
00543   }else{
00544     return 0;
00545   }
00546 }
00547 
00548 void HeavyFlavorValidation::match( MonitorElement * me, vector<LeafCandidate> & from, vector<LeafCandidate> & to, double dRMatchingCut, vector<int> & map ){
00549   vector<double> dR(from.size());
00550   for(size_t i=0; i<from.size(); i++){
00551     map[i] = -1;
00552     dR[i] = 10.;
00553     //find closest
00554     for(size_t j=0; j<to.size(); j++){
00555       double dRtmp = deltaR<double>(from[i].eta(), from[i].phi(), to[j].eta(), to[j].phi());
00556       if( dRtmp < dR[i] ){
00557         dR[i] = dRtmp;
00558         map[i] = j;
00559       }
00560     }
00561     //fill matching histo
00562     if( map[i] != -1 ){
00563       me->Fill( to[map[i]].eta()-from[i].eta(), deltaPhi<double>(to[map[i]].phi(), from[i].phi()) );
00564     }
00565     //apply matching cut
00566     if( dR[i] > dRMatchingCut ){
00567       map[i] = -1;
00568     }
00569     //remove duplication
00570     if( map[i] != -1 ){
00571       for(size_t k=0; k<i; k++){
00572         if( map[k] != -1 && map[i] == map[k] ){
00573           if( dR[i] < dR[k] ){
00574             map[k] = -1;
00575           }else{
00576             map[i] = -1;
00577           }
00578           break;
00579         }
00580       }
00581     }
00582   }
00583 }
00584 
00585 void HeavyFlavorValidation::myBook2D( TString name, vector<double> &ptBins, TString ptLabel, vector<double> &etaBins, TString etaLabel, TString title )
00586 {
00587 //   dqmStore->setCurrentFolder(dqmFolder+"/"+folder);
00588   int ptN = ptBins.size()==3 ? (int)ptBins[0]+1 : ptBins.size();
00589   Double_t *pt = new Double_t[ ptN ];
00590   for(int i=0; i<ptN; i++){
00591     pt[i] = ptBins.size()==3 ? ptBins[1] + i*(ptBins[2]-ptBins[1])/ptBins[0] : ptBins[i] ;
00592   }
00593   int etaN = etaBins.size()==3 ? (int)etaBins[0]+1 : etaBins.size();
00594   Double_t *eta = new Double_t[ etaN ];
00595   for(int i=0; i<etaN; i++){
00596     eta[i] = etaBins.size()==3 ? etaBins[1] + i*(etaBins[2]-etaBins[1])/etaBins[0] : etaBins[i] ;
00597   }
00598   TH2F *h = new TH2F( name, name, ptN-1, pt, etaN-1, eta  );
00599   h->SetXTitle(ptLabel);
00600   h->SetYTitle(etaLabel);
00601   h->SetTitle(title);
00602   ME[name] = dqmStore->book2D( name.Data(), h );
00603   delete h;
00604 }
00605 
00606 void HeavyFlavorValidation::myBookProfile2D( TString name, vector<double> &ptBins, TString ptLabel, vector<double> &etaBins, TString etaLabel, TString title )
00607 {
00608 //   dqmStore->setCurrentFolder(dqmFolder+"/"+folder);
00609   int ptN = ptBins.size()==3 ? (int)ptBins[0]+1 : ptBins.size();
00610   Double_t *pt = new Double_t[ ptN ];
00611   for(int i=0; i<ptN; i++){
00612     pt[i] = ptBins.size()==3 ? ptBins[1] + i*(ptBins[2]-ptBins[1])/ptBins[0] : ptBins[i] ;
00613   }
00614   int etaN = etaBins.size()==3 ? (int)etaBins[0]+1 : etaBins.size();
00615   Double_t *eta = new Double_t[ etaN ];
00616   for(int i=0; i<etaN; i++){
00617     eta[i] = etaBins.size()==3 ? etaBins[1] + i*(etaBins[2]-etaBins[1])/etaBins[0] : etaBins[i] ;
00618   }
00619   TProfile2D *h = new TProfile2D( name, name, ptN-1, pt, etaN-1, eta  );
00620   h->SetXTitle(ptLabel);
00621   h->SetYTitle(etaLabel);
00622   h->SetTitle(title);
00623   ME[name] = dqmStore->bookProfile2D( name.Data(), h );
00624   delete h;
00625 }
00626 
00627 void HeavyFlavorValidation::myBook1D( TString name, vector<double> &bins, TString label, TString title )
00628 {
00629 //   dqmStore->setCurrentFolder(dqmFolder+"/"+folder);
00630   int binsN = bins.size()==3 ? (int)bins[0]+1 : bins.size();
00631   Double_t *myBins = new Double_t[ binsN ];
00632   for(int i=0; i<binsN; i++){
00633     myBins[i] = bins.size()==3 ? bins[1] + i*(bins[2]-bins[1])/bins[0] : bins[i] ;
00634   }
00635   TH1F *h = new TH1F( name, name, binsN-1, myBins );
00636   h->SetXTitle(label);
00637   h->SetTitle(title);
00638   ME[name] = dqmStore->book1D( name.Data(), h );
00639   delete h;
00640 }
00641 
00642 HeavyFlavorValidation::~HeavyFlavorValidation(){
00643 }
00644 
00645 //define this as a plug-in
00646 DEFINE_FWK_MODULE(HeavyFlavorValidation);
/data/refman/pasoursint/CMSSW_6_1_1/src/HLTriggerOffline/HeavyFlavor/src/HeavyFlavorValidation.cc
