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("Filter")){
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           else if(name.Contains("mumuFilter") || name.Contains("JpsiTrackMass"))
00154             level = 4;
00155           filterNamesLevels.push_back( pair<string,int>(moduleNames[j],level) );
00156           os<<" "<<moduleNames[j];
00157         }
00158       }
00159       break;
00160     }
00161   }
00162   if(filterNamesLevels.size()==0){
00163     LogDebug("HLTriggerOfflineHeavyFlavor")<<"Bad Trigger Path: "<<triggerPathName<<endl;
00164     return;
00165   }else{
00166     LogDebug("HLTriggerOfflineHeavyFlavor")<<"Trigger Path: "<<triggerPathName<<" has filters:"<<os.str();
00167   }
00168 
00169 //create Monitor Elements
00170   dqmStore = Service<DQMStore>().operator->();  
00171   if( !dqmStore ){
00172     LogError("HLTriggerOfflineHeavyFlavor") << "Could not find DQMStore service\n";
00173     return;
00174   }
00175   dqmStore->setVerbose(0);
00176   dqmStore->setCurrentFolder((dqmFolder+"/")+triggerProcessName+"/"+triggerPathName);
00177 // Eta Pt Single  
00178   myBook2D( "genMuon_genEtaPt", muonEtaBins, "#mu eta", muonPtBins, " #mu pT (GeV)");
00179   myBook2D( "globMuon_genEtaPt", muonEtaBins, "#mu eta", muonPtBins, " #mu pT (GeV)");
00180   myBook2D( "globMuon_recoEtaPt", muonEtaBins, "#mu eta", muonPtBins, " #mu pT (GeV)");
00181   for(size_t i=0; i<filterNamesLevels.size(); i++){
00182     myBook2D( TString::Format("filt%dMuon_recoEtaPt",int(i+1)), muonEtaBins, "#mu eta", muonPtBins, " #mu pT (GeV)", filterNamesLevels[i].first);
00183   }
00184   myBook2D( "pathMuon_recoEtaPt", muonEtaBins, "#mu eta", muonPtBins, " #mu pT (GeV)", triggerPathName);
00185   myBook2D( "resultMuon_recoEtaPt", muonEtaBins, "#mu eta", muonPtBins, " #mu pT (GeV)");
00186 // Eta Pt Single Resolution
00187   myBookProfile2D( "resGlobGen_genEtaPt", muonEtaBins, "#mu eta", muonPtBins, " #mu pT (GeV)");
00188   for(size_t i=0; i<filterNamesLevels.size(); i++){
00189     myBookProfile2D( TString::Format("resFilt%dGlob_recoEtaPt",int(i+1)), muonEtaBins, "#mu eta", muonPtBins, " #mu pT (GeV)", filterNamesLevels[i].first);
00190   }
00191   myBookProfile2D( "resPathGlob_recoEtaPt", muonEtaBins, "#mu eta", muonPtBins, " #mu pT (GeV)", triggerPathName);
00192 // Eta Pt Double  
00193   myBook2D( "genDimuon_genEtaPt", dimuonEtaBins, "#mu#mu eta", dimuonPtBins, " #mu#mu pT (GeV)");
00194   myBook2D( "globDimuon_genEtaPt", dimuonEtaBins, "#mu#mu eta", dimuonPtBins, " #mu#mu pT (GeV)");
00195   myBook2D( "globDimuon_recoEtaPt", dimuonEtaBins, "#mu#mu eta", dimuonPtBins, " #mu#mu pT (GeV)");
00196   for(size_t i=0; i<filterNamesLevels.size(); i++){
00197     myBook2D( TString::Format("filt%dDimuon_recoEtaPt",int(i+1)), dimuonEtaBins, "#mu#mu eta", dimuonPtBins, " #mu#mu pT (GeV)", filterNamesLevels[i].first);
00198   }  
00199   myBook2D( "pathDimuon_recoEtaPt", dimuonEtaBins, "#mu#mu eta", dimuonPtBins, " #mu#mu pT (GeV)", triggerPathName);
00200   myBook2D( "resultDimuon_recoEtaPt", dimuonEtaBins, "#mu#mu eta", dimuonPtBins, " #mu#mu pT (GeV)");
00201   for(size_t i=0; i<filterNamesLevels.size(); i++){
00202     myBook2D( TString::Format("diFilt%dDimuon_recoEtaPt",int(i+1)), dimuonEtaBins, "#mu#mu eta", dimuonPtBins, " #mu#mu pT (GeV)", filterNamesLevels[i].first);
00203   }  
00204   myBook2D( "diPathDimuon_recoEtaPt", dimuonEtaBins, "#mu#mu eta", dimuonPtBins, " #mu#mu pT (GeV)", triggerPathName);
00205 // Eta Phi Single
00206   myBook2D( "genMuon_genEtaPhi", muonEtaBins, "#mu eta", muonPhiBins, "#mu phi");
00207   myBook2D( "globMuon_genEtaPhi", muonEtaBins, "#mu eta", muonPhiBins, "#mu phi");
00208   myBook2D( "globMuon_recoEtaPhi", muonEtaBins, "#mu eta", muonPhiBins, "#mu phi");
00209   for(size_t i=0; i<filterNamesLevels.size(); i++){
00210     myBook2D( TString::Format("filt%dMuon_recoEtaPhi",int(i+1)), muonEtaBins, "#mu eta", muonPhiBins, "#mu phi", filterNamesLevels[i].first);
00211   }
00212   myBook2D( "pathMuon_recoEtaPhi", muonEtaBins, "#mu eta", muonPhiBins, "#mu phi", triggerPathName);
00213   myBook2D( "resultMuon_recoEtaPhi", muonEtaBins, "#mu eta", muonPhiBins, "#mu phi");
00214 // Rap Pt Double
00215   myBook2D( "genDimuon_genRapPt", dimuonEtaBins, "#mu#mu rapidity", dimuonPtBins, " #mu#mu pT (GeV)");
00216   myBook2D( "globDimuon_genRapPt", dimuonEtaBins, "#mu#mu rapidity", dimuonPtBins, " #mu#mu pT (GeV)");
00217   myBook2D( "globDimuon_recoRapPt", dimuonEtaBins, "#mu#mu rapidity", dimuonPtBins, " #mu#mu pT (GeV)");
00218   for(size_t i=0; i<filterNamesLevels.size(); i++){
00219     myBook2D( TString::Format("filt%dDimuon_recoRapPt",int(i+1)), dimuonEtaBins, "#mu#mu rapidity", dimuonPtBins, " #mu#mu pT (GeV)", filterNamesLevels[i].first);
00220   }  
00221   myBook2D( "pathDimuon_recoRapPt", dimuonEtaBins, "#mu#mu rapidity", dimuonPtBins, " #mu#mu pT (GeV)", triggerPathName);
00222   myBook2D( "resultDimuon_recoRapPt", dimuonEtaBins, "#mu#mu rapidity", dimuonPtBins, " #mu#mu pT (GeV)");
00223   for(size_t i=0; i<filterNamesLevels.size(); i++){
00224     myBook2D( TString::Format("diFilt%dDimuon_recoRapPt",int(i+1)), dimuonEtaBins, "#mu#mu rapidity", dimuonPtBins, " #mu#mu pT (GeV)", filterNamesLevels[i].first);
00225   }  
00226   myBook2D( "diPathDimuon_recoRapPt", dimuonEtaBins, "#mu#mu rapidity", dimuonPtBins, " #mu#mu pT (GeV)", triggerPathName);
00227 // Pt DR Double  
00228   myBook2D( "genDimuon_genPtDR", dimuonPtBins, " #mu#mu pT (GeV)", dimuonDRBins, "#mu#mu #Delta R at IP");
00229   myBook2D( "globDimuon_genPtDR", dimuonPtBins, " #mu#mu pT (GeV)", dimuonDRBins, "#mu#mu #Delta R at IP");
00230   myBook2D( "globDimuon_recoPtDR", dimuonPtBins, " #mu#mu pT (GeV)", dimuonDRBins, "#mu#mu #Delta R at IP");
00231   for(size_t i=0; i<filterNamesLevels.size(); i++){
00232     myBook2D( TString::Format("filt%dDimuon_recoPtDR",int(i+1)), dimuonPtBins, " #mu#mu pT (GeV)", dimuonDRBins, "#mu#mu #Delta R at IP", filterNamesLevels[i].first);
00233   }
00234   myBook2D( "pathDimuon_recoPtDR", dimuonPtBins, " #mu#mu pT (GeV)", dimuonDRBins, "#mu#mu #Delta R at IP", triggerPathName);
00235   for(size_t i=0; i<filterNamesLevels.size(); i++){
00236     myBook2D( TString::Format("diFilt%dDimuon_recoPtDR",int(i+1)), dimuonPtBins, " #mu#mu pT (GeV)", dimuonDRBins, "#mu#mu #Delta R at IP", filterNamesLevels[i].first);
00237   }  
00238   myBook2D( "diPathDimuon_recoPtDR", dimuonPtBins, " #mu#mu pT (GeV)", dimuonDRBins, "#mu#mu #Delta R at IP", triggerPathName);
00239   myBook2D( "resultDimuon_recoPtDR", dimuonPtBins, " #mu#mu pT (GeV)", dimuonDRBins, "#mu#mu #Delta R at IP");
00240 // Pt DRpos Double  
00241 //   myBook2D( "genDimuon_genPtDRpos", dimuonPtBins, " #mu#mu pT (GeV)", dimuonDRBins, "#mu#mu #Delta R in MS");
00242 //   myBook2D( "globDimuon_genPtDRpos", dimuonPtBins, " #mu#mu pT (GeV)", dimuonDRBins, "#mu#mu #Delta R in MS");
00243   myBook2D( "globDimuon_recoPtDRpos", dimuonPtBins, " #mu#mu pT (GeV)", dimuonDRBins, "#mu#mu #Delta R in MS");
00244   for(size_t i=0; i<filterNamesLevels.size(); i++){
00245     myBook2D( TString::Format("filt%dDimuon_recoPtDRpos",int(i+1)), dimuonPtBins, " #mu#mu pT (GeV)", dimuonDRBins, "#mu#mu #Delta R in MS", filterNamesLevels[i].first);
00246   }
00247   myBook2D( "pathDimuon_recoPtDRpos", dimuonPtBins, " #mu#mu pT (GeV)", dimuonDRBins, "#mu#mu #Delta R in MS", triggerPathName);
00248   for(size_t i=0; i<filterNamesLevels.size(); i++){
00249     myBook2D( TString::Format("diFilt%dDimuon_recoPtDRpos",int(i+1)), dimuonPtBins, " #mu#mu pT (GeV)", dimuonDRBins, "#mu#mu #Delta R in MS", filterNamesLevels[i].first);
00250   }  
00251   myBook2D( "diPathDimuon_recoPtDRpos", dimuonPtBins, " #mu#mu pT (GeV)", dimuonDRBins, "#mu#mu #Delta R in MS", triggerPathName);
00252   myBook2D( "resultDimuon_recoPtDRpos", dimuonPtBins, " #mu#mu pT (GeV)", dimuonDRBins, "#mu#mu #Delta R in MS");
00253 
00254 /*  ME["genDimuon_genPt+Pt-"] = myBook2D( "genDimuon_genPt+Pt-", muonPtBins, "#mu+ pt (GeV)", muonPtBins, "#mu- pT (GeV)");
00255   ME["genGlobDimuon_genPt+Pt-"] = myBook2D( "genGlobDimuon_genPt+Pt-", muonPtBins, "#mu+ pt (GeV)", muonPtBins, "#mu- pT (GeV)");
00256   ME["genGlobDimuon_recoPt+Pt-"] = myBook2D( "genGlobDimuon_recoPt+Pt-", muonPtBins, "#mu+ pt (GeV)", muonPtBins, "#mu- pT (GeV)");
00257   ME["genGlobL1Dimuon_recoPt+Pt-"] = myBook2D( "genGlobL1Dimuon_recoPt+Pt-", muonPtBins, "#mu+ pt (GeV)", muonPtBins, "#mu- pT (GeV)");
00258   ME["genGlobL1L2Dimuon_recoPt+Pt-"] = myBook2D( "genGlobL1L2Dimuon_recoPt+Pt-", muonPtBins, "#mu+ pt (GeV)", muonPtBins, "#mu- pT (GeV)");
00259   ME["genGlobL1L2L3Dimuon_recoPt+Pt-"] = myBook2D( "genGlobL1L2L3Dimuon_recoPt+Pt-", muonPtBins, "#mu+ pt (GeV)", muonPtBins, "#mu- pT (GeV)");*/
00260 
00261 // Matching
00262   myBook2D( "globGen_deltaEtaDeltaPhi", deltaEtaBins, "#Delta eta", deltaPhiBins, "#Delta phi");
00263   for(size_t i=0; i<filterNamesLevels.size(); i++){
00264     myBook2D( TString::Format("filt%dGlob_deltaEtaDeltaPhi",int(i+1)), deltaEtaBins, "#Delta eta", deltaPhiBins, "#Delta phi", filterNamesLevels[i].first);
00265   }    
00266   myBook2D( "pathGlob_deltaEtaDeltaPhi", deltaEtaBins, "#Delta eta", deltaPhiBins, "#Delta phi", triggerPathName);
00267 // Size of containers
00268   vector<double> sizeBins; sizeBins.push_back(10); sizeBins.push_back(0); sizeBins.push_back(10);
00269   myBook1D( "genMuon_size", sizeBins, "container size" );
00270   myBook1D( "globMuon_size", sizeBins, "container size" );
00271   for(size_t i=0; i<filterNamesLevels.size(); i++){
00272     myBook1D( TString::Format("filt%dMuon_size",int(i+1)), sizeBins, "container size", filterNamesLevels[i].first);
00273   }    
00274   myBook1D( "pathMuon_size", sizeBins, "container size", triggerPathName );
00275 }
00276 
00277 void HeavyFlavorValidation::analyze(const Event& iEvent, const EventSetup& iSetup){
00278   if( !dqmStore ){
00279     LogDebug("HLTriggerOfflineHeavyFlavor")<<"Could not access DQM Store service"<<endl;
00280     return;
00281   }
00282   if(filterNamesLevels.size()==0){
00283     return;
00284   }
00285 //access the containers and create LeafCandidate copies
00286   vector<LeafCandidate> genMuons;
00287   Handle<GenParticleCollection> genParticles;
00288   iEvent.getByLabel(genParticlesTag, genParticles);
00289   if(genParticles.isValid()){
00290     for(GenParticleCollection::const_iterator p=genParticles->begin(); p!= genParticles->end(); ++p){
00291       if( p->status() == 1 && std::abs(p->pdgId())==13 && 
00292           ( find( motherIDs.begin(), motherIDs.end(), -1 )!=motherIDs.end() || find( motherIDs.begin(), motherIDs.end(), getMotherId(&(*p)) )!=motherIDs.end() ) ){
00293         genMuons.push_back( *p );
00294       }  
00295     }
00296   }else{
00297     LogDebug("HLTriggerOfflineHeavyFlavor")<<"Could not access GenParticleCollection"<<endl;
00298   }
00299   sort(genMuons.begin(), genMuons.end(), GreaterByPt<LeafCandidate>());
00300   ME["genMuon_size"]->Fill(genMuons.size());
00301   LogDebug("HLTriggerOfflineHeavyFlavor")<<"GenParticleCollection from "<<genParticlesTag<<" has size: "<<genMuons.size()<<endl;
00302   
00303   vector<LeafCandidate> globMuons;
00304   vector<LeafCandidate> globMuons_position;
00305   Handle<MuonCollection> recoMuonsHandle;
00306   iEvent.getByLabel(recoMuonsTag, recoMuonsHandle);
00307   if(recoMuonsHandle.isValid()){
00308     for(MuonCollection::const_iterator p=recoMuonsHandle->begin(); p!= recoMuonsHandle->end(); ++p){
00309       if(p->isGlobalMuon()){
00310         globMuons.push_back(*p);
00311         globMuons_position.push_back( LeafCandidate( p->charge(), math::XYZTLorentzVector(p->outerTrack()->innerPosition().x(), p->outerTrack()->innerPosition().y(), p->outerTrack()->innerPosition().z(), 0.) ) );
00312       }
00313     }
00314   }else{
00315     LogDebug("HLTriggerOfflineHeavyFlavor")<<"Could not access reco Muons"<<endl;
00316   }
00317   ME["globMuon_size"]->Fill(globMuons.size());
00318   LogDebug("HLTriggerOfflineHeavyFlavor")<<"Global Muons from "<<recoMuonsTag<<" has size: "<<globMuons.size()<<endl;
00319 
00320 // access RAW trigger event  
00321   vector<vector<LeafCandidate> > muonsAtFilter;
00322   vector<vector<LeafCandidate> > muonPositionsAtFilter;  
00323   for(size_t i=0; i<filterNamesLevels.size(); i++){
00324     muonsAtFilter.push_back(vector<LeafCandidate>());
00325     muonPositionsAtFilter.push_back(vector<LeafCandidate>());
00326   }
00327   Handle<TriggerEventWithRefs> rawTriggerEvent;
00328   iEvent.getByLabel( triggerSummaryRAWTag, rawTriggerEvent );
00329   if( rawTriggerEvent.isValid() ){
00330     for(size_t i=0; i<filterNamesLevels.size(); i++){
00331       size_t index = rawTriggerEvent->filterIndex(InputTag( filterNamesLevels[i].first, "", triggerProcessName ));
00332       if ( index < rawTriggerEvent->size() ){
00333         if( filterNamesLevels[i].second==1 ){
00334           vector<L1MuonParticleRef> l1Cands;
00335           rawTriggerEvent->getObjects( index, TriggerL1Mu, l1Cands );
00336           for(size_t j=0; j<l1Cands.size(); j++){
00337             muonsAtFilter[i].push_back(*l1Cands[j]);
00338           }
00339         }else{
00340           vector<RecoChargedCandidateRef> hltCands;        
00341           rawTriggerEvent->getObjects( index, TriggerMuon, hltCands );
00342           for(size_t j=0; j<hltCands.size(); j++){
00343             muonsAtFilter[i].push_back(*hltCands[j]);
00344             if( filterNamesLevels[i].second==2 ){
00345               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.) ) );
00346             }
00347           }
00348         }
00349       }
00350       ME[TString::Format("filt%dMuon_size",int(i+1))]->Fill(muonsAtFilter[i].size());
00351       LogDebug("HLTriggerOfflineHeavyFlavor")<<"Filter \""<<filterNamesLevels[i].first<<"\" has "<<muonsAtFilter[i].size()<<" muons"<<endl;
00352     }
00353   }else{
00354     LogDebug("HLTriggerOfflineHeavyFlavor")<<"Could not access RAWTriggerEvent"<<endl;
00355   }
00356 
00357 // access AOD trigger event  
00358   vector<LeafCandidate> pathMuons;
00359   Handle<TriggerEvent> aodTriggerEvent;
00360   iEvent.getByLabel(triggerSummaryAODTag,aodTriggerEvent);
00361   if(aodTriggerEvent.isValid()){
00362     TriggerObjectCollection allObjects = aodTriggerEvent->getObjects();
00363     for(int i=0; i<aodTriggerEvent->sizeFilters(); i++){         
00364      if(aodTriggerEvent->filterTag(i)==InputTag((filterNamesLevels.end()-1)->first,"",triggerProcessName)){
00365         Keys keys = aodTriggerEvent->filterKeys(i);
00366         for(size_t j=0; j<keys.size(); j++){
00367           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 ) ) );
00368         }
00369       }
00370     }
00371     ME["pathMuon_size"]->Fill(pathMuons.size());
00372     LogDebug("HLTriggerOfflineHeavyFlavor")<<"Path \""<<triggerPathName<<"\" has "<<pathMuons.size()<<" muons at last filter \""<<(filterNamesLevels.end()-1)->first<<"\""<<endl;
00373   }else{
00374     LogDebug("HLTriggerOfflineHeavyFlavor")<<"Could not access AODTriggerEvent"<<endl;
00375   }
00376 
00377 // access Trigger Results
00378   bool triggerFired = false;
00379   Handle<TriggerResults> triggerResults;
00380   iEvent.getByLabel(triggerResultsTag,triggerResults);
00381   if(triggerResults.isValid()){
00382     LogDebug("HLTriggerOfflineHeavyFlavor")<<"Successfully initialized "<<triggerResultsTag<<endl;
00383     const edm::TriggerNames & triggerNames = iEvent.triggerNames(*triggerResults);
00384     bool hlt_exists = false;
00385     for ( unsigned int i=0; i!=triggerNames.size(); i++) {
00386       TString hlt_name = triggerNames.triggerName(i);
00387       if (hlt_name.Contains(triggerPathName)) {
00388         triggerFired = triggerResults->accept( i );
00389         hlt_exists = true;
00390         break;
00391       }
00392     }
00393     if (!hlt_exists) {
00394       LogDebug("HLTriggerOfflineHeavyFlavor")<<triggerResultsTag<<" has no trigger: "<<triggerPathName<<endl;
00395     }
00396   }else{
00397     LogDebug("HLTriggerOfflineHeavyFlavor")<<"Could not initialize "<<triggerResultsTag<<endl;
00398   }
00399 
00400 //create matching maps
00401   vector<int> glob_gen(genMuons.size(),-1);
00402   match( ME["globGen_deltaEtaDeltaPhi"], genMuons, globMuons, genGlobDeltaRMatchingCut, glob_gen );
00403   vector<vector<int> > filt_glob;
00404   for(size_t i=0; i<filterNamesLevels.size(); i++){
00405     filt_glob.push_back( vector<int>(globMuons.size(),-1) );            
00406     if( filterNamesLevels[i].second == 1 ){
00407       match( ME[TString::Format("filt%dGlob_deltaEtaDeltaPhi",int(i+1))], globMuons_position, muonsAtFilter[i] ,globL1DeltaRMatchingCut, filt_glob[i] );
00408     }else if( filterNamesLevels[i].second == 2 ){
00409       match( ME[TString::Format("filt%dGlob_deltaEtaDeltaPhi",int(i+1))], globMuons_position, muonPositionsAtFilter[i] ,globL2DeltaRMatchingCut, filt_glob[i] );
00410     }else if( filterNamesLevels[i].second == 3 || filterNamesLevels[i].second == 4){
00411       match( ME[TString::Format("filt%dGlob_deltaEtaDeltaPhi",int(i+1))], globMuons, muonsAtFilter[i] ,globL3DeltaRMatchingCut, filt_glob[i] );
00412     }
00413   }
00414   vector<int> path_glob(globMuons.size(),-1);
00415   if( (filterNamesLevels.end()-1)->second == 1 ){
00416     match( ME["pathGlob_deltaEtaDeltaPhi"], globMuons_position, pathMuons ,globL1DeltaRMatchingCut, path_glob );
00417   }else if( (filterNamesLevels.end()-1)->second == 2 ){
00418     match( ME["pathGlob_deltaEtaDeltaPhi"], globMuons, pathMuons ,globL2DeltaRMatchingCut, path_glob );
00419   }else if( (filterNamesLevels.end()-1)->second == 3 || (filterNamesLevels.end()-1)->second == 4){
00420     match( ME["pathGlob_deltaEtaDeltaPhi"], globMuons, pathMuons ,globL3DeltaRMatchingCut, path_glob );
00421   }
00422     
00423 //fill histos
00424   bool first = true;
00425   for(size_t i=0; i<genMuons.size(); i++){
00426     ME["genMuon_genEtaPt"]->Fill(genMuons[i].eta(), genMuons[i].pt());
00427     ME["genMuon_genEtaPhi"]->Fill(genMuons[i].eta(), genMuons[i].phi());
00428     if(glob_gen[i] != -1){
00429       ME["resGlobGen_genEtaPt"]->Fill(genMuons[i].eta(), genMuons[i].pt(), (globMuons[glob_gen[i]].pt()-genMuons[i].pt())/genMuons[i].pt() );
00430       ME["globMuon_genEtaPt"]->Fill(genMuons[i].eta(), genMuons[i].pt());
00431       ME["globMuon_genEtaPhi"]->Fill(genMuons[i].eta(), genMuons[i].phi());
00432       ME["globMuon_recoEtaPt"]->Fill(globMuons[glob_gen[i]].eta(), globMuons[glob_gen[i]].pt());
00433       ME["globMuon_recoEtaPhi"]->Fill(globMuons[glob_gen[i]].eta(), globMuons[glob_gen[i]].phi());
00434       for(size_t f=0; f<filterNamesLevels.size(); f++){
00435         if(filt_glob[f][glob_gen[i]] != -1){
00436           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() );
00437           ME[TString::Format("filt%dMuon_recoEtaPt",int(f+1))]->Fill(globMuons[glob_gen[i]].eta(), globMuons[glob_gen[i]].pt());
00438           ME[TString::Format("filt%dMuon_recoEtaPhi",int(f+1))]->Fill(globMuons[glob_gen[i]].eta(), globMuons[glob_gen[i]].phi());
00439         }else{
00440           break;
00441         }  
00442       }
00443       if(path_glob[glob_gen[i]] != -1){
00444         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() );
00445         ME["pathMuon_recoEtaPt"]->Fill(globMuons[glob_gen[i]].eta(), globMuons[glob_gen[i]].pt());
00446         ME["pathMuon_recoEtaPhi"]->Fill(globMuons[glob_gen[i]].eta(), globMuons[glob_gen[i]].phi());
00447       }
00448 //highest pt muon      
00449       if( first ){
00450         first = false;
00451         if( triggerFired ){
00452           ME["resultMuon_recoEtaPt"]->Fill(globMuons[glob_gen[i]].eta(), globMuons[glob_gen[i]].pt());
00453           ME["resultMuon_recoEtaPhi"]->Fill(globMuons[glob_gen[i]].eta(), globMuons[glob_gen[i]].phi());
00454         }
00455       }
00456     }
00457   }  
00458 
00459 //fill dimuon histograms (highest pT, opposite charge) 
00460   int secondMuon = 0;
00461   for(size_t j=1; j<genMuons.size(); j++){
00462     if(genMuons[0].charge()*genMuons[j].charge()==-1){
00463       secondMuon = j;
00464       break;
00465     }
00466   }
00467   if(secondMuon > 0){
00468 //    int pos = genMuons[0].charge()>0 ? 0 : secondMuon ;
00469 //    int neg = genMuons[0].charge()<0 ? 0 : secondMuon ;
00470 //two generated
00471     double genDimuonPt = (genMuons[0].p4()+genMuons[secondMuon].p4()).pt();
00472     double genDimuonEta = (genMuons[0].p4()+genMuons[secondMuon].p4()).eta();
00473     double genDimuonRap = (genMuons[0].p4()+genMuons[secondMuon].p4()).Rapidity();
00474     double genDimuonDR = deltaR<LeafCandidate,LeafCandidate>( genMuons[0], genMuons[secondMuon] );
00475     bool highPt = genMuons[0].pt()>7. && genMuons[secondMuon].pt()>7;
00476     ME["genDimuon_genEtaPt"]->Fill( genDimuonEta, genDimuonPt );
00477     ME["genDimuon_genRapPt"]->Fill( genDimuonRap, genDimuonPt );
00478     if(highPt) ME["genDimuon_genPtDR"]->Fill( genDimuonPt, genDimuonDR );
00479 //two global
00480     if(glob_gen[0]!=-1 && glob_gen[secondMuon]!=-1){
00481       ME["globDimuon_genEtaPt"]->Fill( genDimuonEta, genDimuonPt );
00482       ME["globDimuon_genRapPt"]->Fill( genDimuonRap, genDimuonPt );
00483       if(highPt) ME["globDimuon_genPtDR"]->Fill( genDimuonPt, genDimuonDR );
00484       double globDimuonPt = (globMuons[glob_gen[0]].p4()+globMuons[glob_gen[secondMuon]].p4()).pt();
00485       double globDimuonEta = (globMuons[glob_gen[0]].p4()+globMuons[glob_gen[secondMuon]].p4()).eta();
00486       double globDimuonRap = (globMuons[glob_gen[0]].p4()+globMuons[glob_gen[secondMuon]].p4()).Rapidity();
00487       double globDimuonDR = deltaR<LeafCandidate,LeafCandidate>( globMuons[glob_gen[0]], globMuons[glob_gen[secondMuon]] );
00488       double globDimuonDRpos = deltaR<LeafCandidate,LeafCandidate>( globMuons_position[glob_gen[0]], globMuons_position[glob_gen[secondMuon]] );
00489       ME["globDimuon_recoEtaPt"]->Fill( globDimuonEta, globDimuonPt );
00490       ME["globDimuon_recoRapPt"]->Fill( globDimuonRap, globDimuonPt );
00491       if(highPt) ME["globDimuon_recoPtDR"]->Fill( globDimuonPt, globDimuonDR );
00492       if(highPt) ME["globDimuon_recoPtDRpos"]->Fill( globDimuonPt, globDimuonDRpos );
00493 //two filter objects    
00494       for(size_t f=0; f<filterNamesLevels.size(); f++){
00495         if(filt_glob[f][glob_gen[0]] != -1 && filt_glob[f][glob_gen[secondMuon]] != -1){
00496           ME[TString::Format("diFilt%dDimuon_recoEtaPt",int(f+1))]->Fill( globDimuonEta, globDimuonPt );
00497           ME[TString::Format("diFilt%dDimuon_recoRapPt",int(f+1))]->Fill( globDimuonRap, globDimuonPt );
00498           if(highPt) ME[TString::Format("diFilt%dDimuon_recoPtDR",int(f+1))]->Fill( globDimuonPt, globDimuonDR );
00499           if(highPt) ME[TString::Format("diFilt%dDimuon_recoPtDRpos",int(f+1))]->Fill( globDimuonPt, globDimuonDRpos );
00500         }else{
00501           break;
00502         }  
00503       } 
00504 //one filter object
00505       for(size_t f=0; f<filterNamesLevels.size(); f++){
00506         if(filt_glob[f][glob_gen[0]] != -1 || filt_glob[f][glob_gen[secondMuon]] != -1){
00507           ME[TString::Format("filt%dDimuon_recoEtaPt",int(f+1))]->Fill( globDimuonEta, globDimuonPt );
00508           ME[TString::Format("filt%dDimuon_recoRapPt",int(f+1))]->Fill( globDimuonRap, globDimuonPt );
00509           if(highPt) ME[TString::Format("filt%dDimuon_recoPtDR",int(f+1))]->Fill( globDimuonPt, globDimuonDR );
00510           if(highPt) ME[TString::Format("filt%dDimuon_recoPtDRpos",int(f+1))]->Fill( globDimuonPt, globDimuonDRpos );
00511         }else{
00512           break;
00513         }
00514       }
00515 //two path objects
00516       if(path_glob[glob_gen[0]] != -1 && path_glob[glob_gen[secondMuon]] != -1){
00517         ME["diPathDimuon_recoEtaPt"]->Fill( globDimuonEta, globDimuonPt );
00518         ME["diPathDimuon_recoRapPt"]->Fill( globDimuonRap, globDimuonPt );
00519         if(highPt) ME["diPathDimuon_recoPtDR"]->Fill( globDimuonPt, globDimuonDR );
00520         if(highPt) ME["diPathDimuon_recoPtDRpos"]->Fill( globDimuonPt, globDimuonDRpos );
00521       }
00522 //one path object
00523       if(path_glob[glob_gen[0]] != -1 || path_glob[glob_gen[secondMuon]] != -1){
00524         ME["pathDimuon_recoEtaPt"]->Fill( globDimuonEta, globDimuonPt );
00525         ME["pathDimuon_recoRapPt"]->Fill( globDimuonRap, globDimuonPt );
00526         if(highPt) ME["pathDimuon_recoPtDR"]->Fill( globDimuonPt, globDimuonDR );
00527         if(highPt) ME["pathDimuon_recoPtDRpos"]->Fill( globDimuonPt, globDimuonDRpos );
00528       }
00529 //trigger result
00530       if( triggerFired ){
00531         ME["resultDimuon_recoEtaPt"]->Fill( globDimuonEta, globDimuonPt );
00532         ME["resultDimuon_recoRapPt"]->Fill( globDimuonRap, globDimuonPt );
00533         if(highPt) ME["resultDimuon_recoPtDR"]->Fill( globDimuonPt, globDimuonDR );
00534         if(highPt) ME["resultDimuon_recoPtDRpos"]->Fill( globDimuonPt, globDimuonDRpos );
00535       }
00536     }
00537   }
00538 }
00539 
00540 void HeavyFlavorValidation::endJob(){
00541 }
00542 
00543 int HeavyFlavorValidation::getMotherId( const Candidate * p ){
00544   const Candidate* mother = p->mother();
00545   if( mother ){
00546     if( mother->pdgId() == p->pdgId() ){
00547       return getMotherId(mother);
00548     }else{
00549       return mother->pdgId();
00550     }
00551   }else{
00552     return 0;
00553   }
00554 }
00555 
00556 void HeavyFlavorValidation::match( MonitorElement * me, vector<LeafCandidate> & from, vector<LeafCandidate> & to, double dRMatchingCut, vector<int> & map ){
00557   vector<double> dR(from.size());
00558   for(size_t i=0; i<from.size(); i++){
00559     map[i] = -1;
00560     dR[i] = 10.;
00561     //find closest
00562     for(size_t j=0; j<to.size(); j++){
00563       double dRtmp = deltaR<double>(from[i].eta(), from[i].phi(), to[j].eta(), to[j].phi());
00564       if( dRtmp < dR[i] ){
00565         dR[i] = dRtmp;
00566         map[i] = j;
00567       }
00568     }
00569     //fill matching histo
00570     if( map[i] != -1 ){
00571       me->Fill( to[map[i]].eta()-from[i].eta(), deltaPhi<double>(to[map[i]].phi(), from[i].phi()) );
00572     }
00573     //apply matching cut
00574     if( dR[i] > dRMatchingCut ){
00575       map[i] = -1;
00576     }
00577     //remove duplication
00578     if( map[i] != -1 ){
00579       for(size_t k=0; k<i; k++){
00580         if( map[k] != -1 && map[i] == map[k] ){
00581           if( dR[i] < dR[k] ){
00582             map[k] = -1;
00583           }else{
00584             map[i] = -1;
00585           }
00586           break;
00587         }
00588       }
00589     }
00590   }
00591 }
00592 
00593 void HeavyFlavorValidation::myBook2D( TString name, vector<double> &ptBins, TString ptLabel, vector<double> &etaBins, TString etaLabel, TString title )
00594 {
00595 //   dqmStore->setCurrentFolder(dqmFolder+"/"+folder);
00596   int ptN = ptBins.size()==3 ? (int)ptBins[0]+1 : ptBins.size();
00597   Double_t *pt = new Double_t[ ptN ];
00598   for(int i=0; i<ptN; i++){
00599     pt[i] = ptBins.size()==3 ? ptBins[1] + i*(ptBins[2]-ptBins[1])/ptBins[0] : ptBins[i] ;
00600   }
00601   int etaN = etaBins.size()==3 ? (int)etaBins[0]+1 : etaBins.size();
00602   Double_t *eta = new Double_t[ etaN ];
00603   for(int i=0; i<etaN; i++){
00604     eta[i] = etaBins.size()==3 ? etaBins[1] + i*(etaBins[2]-etaBins[1])/etaBins[0] : etaBins[i] ;
00605   }
00606   TH2F *h = new TH2F( name, name, ptN-1, pt, etaN-1, eta  );
00607   h->SetXTitle(ptLabel);
00608   h->SetYTitle(etaLabel);
00609   h->SetTitle(title);
00610   ME[name] = dqmStore->book2D( name.Data(), h );
00611   delete h;
00612 }
00613 
00614 void HeavyFlavorValidation::myBookProfile2D( TString name, vector<double> &ptBins, TString ptLabel, vector<double> &etaBins, TString etaLabel, TString title )
00615 {
00616 //   dqmStore->setCurrentFolder(dqmFolder+"/"+folder);
00617   int ptN = ptBins.size()==3 ? (int)ptBins[0]+1 : ptBins.size();
00618   Double_t *pt = new Double_t[ ptN ];
00619   for(int i=0; i<ptN; i++){
00620     pt[i] = ptBins.size()==3 ? ptBins[1] + i*(ptBins[2]-ptBins[1])/ptBins[0] : ptBins[i] ;
00621   }
00622   int etaN = etaBins.size()==3 ? (int)etaBins[0]+1 : etaBins.size();
00623   Double_t *eta = new Double_t[ etaN ];
00624   for(int i=0; i<etaN; i++){
00625     eta[i] = etaBins.size()==3 ? etaBins[1] + i*(etaBins[2]-etaBins[1])/etaBins[0] : etaBins[i] ;
00626   }
00627   TProfile2D *h = new TProfile2D( name, name, ptN-1, pt, etaN-1, eta  );
00628   h->SetXTitle(ptLabel);
00629   h->SetYTitle(etaLabel);
00630   h->SetTitle(title);
00631   ME[name] = dqmStore->bookProfile2D( name.Data(), h );
00632   delete h;
00633 }
00634 
00635 void HeavyFlavorValidation::myBook1D( TString name, vector<double> &bins, TString label, TString title )
00636 {
00637 //   dqmStore->setCurrentFolder(dqmFolder+"/"+folder);
00638   int binsN = bins.size()==3 ? (int)bins[0]+1 : bins.size();
00639   Double_t *myBins = new Double_t[ binsN ];
00640   for(int i=0; i<binsN; i++){
00641     myBins[i] = bins.size()==3 ? bins[1] + i*(bins[2]-bins[1])/bins[0] : bins[i] ;
00642   }
00643   TH1F *h = new TH1F( name, name, binsN-1, myBins );
00644   h->SetXTitle(label);
00645   h->SetTitle(title);
00646   ME[name] = dqmStore->book1D( name.Data(), h );
00647   delete h;
00648 }
00649 
00650 HeavyFlavorValidation::~HeavyFlavorValidation(){
00651 }
00652 
00653 //define this as a plug-in
00654 DEFINE_FWK_MODULE(HeavyFlavorValidation);
/afs/cern.ch/work/a/aaltunda/public/www/CMSSW_6_2_5/src/HLTriggerOffline/HeavyFlavor/src/HeavyFlavorValidation.cc
