CMS 3D CMS Logo

 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
HeavyFlavorValidation.cc
Go to the documentation of this file.
1 
8 // Original Author: Zoltan Gecse
9 
10 #include <memory>
20 
36 
38 
42 
44 
45 #include "TLorentzVector.h"
46 
47 using namespace std;
48 using namespace edm;
49 using namespace reco;
50 using namespace l1extra;
51 using namespace trigger;
52 
54  public:
57  protected:
58  void dqmBeginRun(const edm::Run&, const edm::EventSetup&);
59  void bookHistograms(DQMStore::IBooker &, edm::Run const &, edm::EventSetup const &) override;
60  virtual void analyze(const edm::Event&, const edm::EventSetup&) override;
61  private:
62  int getMotherId( const Candidate * p );
63  void match( MonitorElement * me, vector<LeafCandidate> & from, vector<LeafCandidate> & to, double deltaRMatchingCut, vector<int> & map );
64  void myBook2D(DQMStore::IBooker & ibooker, TString name, vector<double> &xBins, TString xLabel, vector<double> &yBins, TString yLabel, TString title);
65  void myBook2D(DQMStore::IBooker & ibooker, TString name, vector<double> &xBins, TString xLabel, vector<double> &yBins, TString yLabel){
66  myBook2D(ibooker, name, xBins, xLabel, yBins, yLabel, name);
67  }
68  void myBookProfile2D(DQMStore::IBooker & ibooker, TString name, vector<double> &xBins, TString xLabel, vector<double> &yBins, TString yLabel, TString title);
69  void myBookProfile2D(DQMStore::IBooker & ibooker, TString name, vector<double> &xBins, TString xLabel, vector<double> &yBins, TString yLabel){
70  myBookProfile2D(ibooker, name, xBins, xLabel, yBins, yLabel, name);
71  }
72  void myBook1D(DQMStore::IBooker & ibooker, TString name, vector<double> &xBins, TString label, TString title );
73  void myBook1D(DQMStore::IBooker & ibooker, TString name, vector<double> &xBins, TString label ){
74  myBook1D(ibooker, name, xBins, label, name );
75  }
76  string dqmFolder;
79 
88 
89  vector<int> motherIDs;
94  vector<double> deltaEtaBins;
95  vector<double> deltaPhiBins;
96  vector<double> muonPtBins;
97  vector<double> muonEtaBins;
98  vector<double> muonPhiBins;
99  vector<double> dimuonPtBins;
100  vector<double> dimuonEtaBins;
101  vector<double> dimuonDRBins;
103  map<TString, MonitorElement *> ME;
104  vector<pair<string,int> > filterNamesLevels;
105  const double muonMass;
106 };
107 
109 //get parameters
110  dqmFolder(pset.getUntrackedParameter<string>("DQMFolder")),
111  triggerProcessName(pset.getUntrackedParameter<string>("TriggerProcessName")),
112  triggerPathName(pset.getUntrackedParameter<string>("TriggerPathName")),
113  motherIDs(pset.getUntrackedParameter<vector<int> >("MotherIDs")),
114  genGlobDeltaRMatchingCut(pset.getUntrackedParameter<double>("GenGlobDeltaRMatchingCut")),
115  globL1DeltaRMatchingCut(pset.getUntrackedParameter<double>("GlobL1DeltaRMatchingCut")),
116  globL2DeltaRMatchingCut(pset.getUntrackedParameter<double>("GlobL2DeltaRMatchingCut")),
117  globL3DeltaRMatchingCut(pset.getUntrackedParameter<double>("GlobL3DeltaRMatchingCut")),
118  deltaEtaBins(pset.getUntrackedParameter<vector<double> >("DeltaEtaBins")),
119  deltaPhiBins(pset.getUntrackedParameter<vector<double> >("DeltaPhiBins")),
120  muonPtBins(pset.getUntrackedParameter<vector<double> >("MuonPtBins")),
121  muonEtaBins(pset.getUntrackedParameter<vector<double> >("MuonEtaBins")),
122  muonPhiBins(pset.getUntrackedParameter<vector<double> >("MuonPhiBins")),
123  dimuonPtBins(pset.getUntrackedParameter<vector<double> >("DimuonPtBins")),
124  dimuonEtaBins(pset.getUntrackedParameter<vector<double> >("DimuonEtaBins")),
125  dimuonDRBins(pset.getUntrackedParameter<vector<double> >("DimuonDRBins")),
126  muonMass(0.106)
127 {
128  triggerSummaryRAWTag = consumes<TriggerEventWithRefs>(InputTag( pset.getUntrackedParameter<string>("TriggerSummaryRAW"), "", triggerProcessName));
129  triggerSummaryAODTag = consumes<TriggerEvent>(InputTag( pset.getUntrackedParameter<string>("TriggerSummaryAOD"), "", triggerProcessName));
130  triggerResultsTag = InputTag( pset.getUntrackedParameter<string>("TriggerResults"), "", triggerProcessName);
131  triggerResultsToken = consumes<TriggerResults>(triggerResultsTag);
132  recoMuonsTag = pset.getParameter<InputTag>("RecoMuons");
133  recoMuonsToken = consumes<MuonCollection>(recoMuonsTag);
134  genParticlesTag = pset.getParameter<InputTag>("GenParticles");
135  genParticlesToken = consumes<GenParticleCollection>(genParticlesTag);
136 }
137 
138 void HeavyFlavorValidation::dqmBeginRun(const edm::Run& iRun, const edm::EventSetup & iSetup ) {
139  //discover HLT configuration
141  bool isChanged;
142  if(hltConfig.init(iRun, iSetup, triggerProcessName, isChanged)){
143  LogDebug("HLTriggerOfflineHeavyFlavor") << "Successfully initialized HLTConfigProvider with process name: "<<triggerProcessName<<endl;
144  }else{
145  LogWarning("HLTriggerOfflineHeavyFlavor") << "Could not initialize HLTConfigProvider with process name: "<<triggerProcessName<<endl;
146  return;
147  }
148  stringstream os;
149  vector<string> triggerNames = hltConfig.triggerNames();
150  for( size_t i = 0; i < triggerNames.size(); i++) {
151  TString triggerName = triggerNames[i];
152  if (triggerName.Contains(triggerPathName)){
153  vector<string> moduleNames = hltConfig.moduleLabels( triggerNames[i] );
154  for( size_t j = 0; j < moduleNames.size(); j++) {
155  TString name = moduleNames[j];
156  if(name.Contains("Filter")){
157  int level = 0;
158  if(name.Contains("L1"))
159  level = 1;
160  if(name.Contains("L2"))
161  level = 2;
162  if(name.Contains("L3"))
163  level = 3;
164  if(name.Contains("mumuFilter") || name.Contains("DiMuon") || name.Contains("MuonL3Filtered") || name.Contains("TrackMassFiltered"))
165  level = 4;
166  if(name.Contains("Vertex") || name.Contains("Dz"))
167  level = 5;
168  filterNamesLevels.push_back( pair<string,int>(moduleNames[j],level) );
169  os<<" "<<moduleNames[j];
170  }
171  }
172  break;
173  }
174  }
175  if(filterNamesLevels.size()==0){
176  LogDebug("HLTriggerOfflineHeavyFlavor")<<"Bad Trigger Path: "<<triggerPathName<<endl;
177  return;
178  }else{
179  LogDebug("HLTriggerOfflineHeavyFlavor")<<"Trigger Path: "<<triggerPathName<<" has filters:"<<os.str();
180  }
181 }
182 
183 
184 
186  edm::Run const & iRun,
187  edm::EventSetup const & iSetup) {
188  ibooker.cd();
190 
191  // create Monitor Elements
192  // Eta Pt Single
193  myBook2D(ibooker, "genMuon_genEtaPt", muonEtaBins, "#mu eta", muonPtBins, " #mu pT (GeV)");
194  myBook2D(ibooker, "globMuon_genEtaPt", muonEtaBins, "#mu eta", muonPtBins, " #mu pT (GeV)");
195  myBook2D(ibooker, "globMuon_recoEtaPt", muonEtaBins, "#mu eta", muonPtBins, " #mu pT (GeV)");
196 
197  for(size_t i=0; i<filterNamesLevels.size(); i++){
198  myBook2D(ibooker, TString::Format("filt%dMuon_recoEtaPt",int(i+1)), muonEtaBins, "#mu eta", muonPtBins, " #mu pT (GeV)", filterNamesLevels[i].first);
199  }
200  myBook2D(ibooker, "pathMuon_recoEtaPt", muonEtaBins, "#mu eta", muonPtBins, " #mu pT (GeV)", triggerPathName);
201  myBook2D(ibooker, "resultMuon_recoEtaPt", muonEtaBins, "#mu eta", muonPtBins, " #mu pT (GeV)");
202  // Eta Pt Single Resolution
203  myBookProfile2D(ibooker, "resGlobGen_genEtaPt", muonEtaBins, "#mu eta", muonPtBins, " #mu pT (GeV)");
204  for(size_t i=0; i<filterNamesLevels.size(); i++){
205  myBookProfile2D(ibooker, TString::Format("resFilt%dGlob_recoEtaPt",int(i+1)), muonEtaBins, "#mu eta", muonPtBins, " #mu pT (GeV)", filterNamesLevels[i].first);
206  }
207  myBookProfile2D(ibooker, "resPathGlob_recoEtaPt", muonEtaBins, "#mu eta", muonPtBins, " #mu pT (GeV)", triggerPathName);
208  // Eta Pt Double
209  myBook2D(ibooker, "genDimuon_genEtaPt", dimuonEtaBins, "#mu#mu eta", dimuonPtBins, " #mu#mu pT (GeV)");
210  myBook2D(ibooker, "globDimuon_genEtaPt", dimuonEtaBins, "#mu#mu eta", dimuonPtBins, " #mu#mu pT (GeV)");
211  myBook2D(ibooker, "globDimuon_recoEtaPt", dimuonEtaBins, "#mu#mu eta", dimuonPtBins, " #mu#mu pT (GeV)");
212  for(size_t i=0; i<filterNamesLevels.size(); i++){
213  myBook2D(ibooker, TString::Format("filt%dDimuon_recoEtaPt",int(i+1)), dimuonEtaBins, "#mu#mu eta", dimuonPtBins, " #mu#mu pT (GeV)", filterNamesLevels[i].first);
214  }
215  myBook2D(ibooker, "pathDimuon_recoEtaPt", dimuonEtaBins, "#mu#mu eta", dimuonPtBins, " #mu#mu pT (GeV)", triggerPathName);
216  myBook2D(ibooker, "resultDimuon_recoEtaPt", dimuonEtaBins, "#mu#mu eta", dimuonPtBins, " #mu#mu pT (GeV)");
217  for(size_t i=0; i<filterNamesLevels.size(); i++){
218  myBook2D(ibooker, TString::Format("diFilt%dDimuon_recoEtaPt",int(i+1)), dimuonEtaBins, "#mu#mu eta", dimuonPtBins, " #mu#mu pT (GeV)", filterNamesLevels[i].first);
219  }
220  myBook2D(ibooker, "diPathDimuon_recoEtaPt", dimuonEtaBins, "#mu#mu eta", dimuonPtBins, " #mu#mu pT (GeV)", triggerPathName);
221  // Eta Phi Single
222  myBook2D(ibooker, "genMuon_genEtaPhi", muonEtaBins, "#mu eta", muonPhiBins, "#mu phi");
223  myBook2D(ibooker, "globMuon_genEtaPhi", muonEtaBins, "#mu eta", muonPhiBins, "#mu phi");
224  myBook2D(ibooker, "globMuon_recoEtaPhi", muonEtaBins, "#mu eta", muonPhiBins, "#mu phi");
225  for(size_t i=0; i<filterNamesLevels.size(); i++){
226  myBook2D(ibooker, TString::Format("filt%dMuon_recoEtaPhi",int(i+1)), muonEtaBins, "#mu eta", muonPhiBins, "#mu phi", filterNamesLevels[i].first);
227  }
228  myBook2D(ibooker, "pathMuon_recoEtaPhi", muonEtaBins, "#mu eta", muonPhiBins, "#mu phi", triggerPathName);
229  myBook2D(ibooker, "resultMuon_recoEtaPhi", muonEtaBins, "#mu eta", muonPhiBins, "#mu phi");
230  // Rap Pt Double
231  myBook2D(ibooker, "genDimuon_genRapPt", dimuonEtaBins, "#mu#mu rapidity", dimuonPtBins, " #mu#mu pT (GeV)");
232  myBook2D(ibooker, "globDimuon_genRapPt", dimuonEtaBins, "#mu#mu rapidity", dimuonPtBins, " #mu#mu pT (GeV)");
233  myBook2D(ibooker, "globDimuon_recoRapPt", dimuonEtaBins, "#mu#mu rapidity", dimuonPtBins, " #mu#mu pT (GeV)");
234  for(size_t i=0; i<filterNamesLevels.size(); i++){
235  myBook2D(ibooker, TString::Format("filt%dDimuon_recoRapPt",int(i+1)), dimuonEtaBins, "#mu#mu rapidity", dimuonPtBins, " #mu#mu pT (GeV)", filterNamesLevels[i].first);
236  }
237  myBook2D(ibooker, "pathDimuon_recoRapPt", dimuonEtaBins, "#mu#mu rapidity", dimuonPtBins, " #mu#mu pT (GeV)", triggerPathName);
238  myBook2D(ibooker, "resultDimuon_recoRapPt", dimuonEtaBins, "#mu#mu rapidity", dimuonPtBins, " #mu#mu pT (GeV)");
239  for(size_t i=0; i<filterNamesLevels.size(); i++){
240  myBook2D(ibooker, TString::Format("diFilt%dDimuon_recoRapPt",int(i+1)), dimuonEtaBins, "#mu#mu rapidity", dimuonPtBins, " #mu#mu pT (GeV)", filterNamesLevels[i].first);
241  }
242  myBook2D(ibooker, "diPathDimuon_recoRapPt", dimuonEtaBins, "#mu#mu rapidity", dimuonPtBins, " #mu#mu pT (GeV)", triggerPathName);
243  // Pt DR Double
244  myBook2D(ibooker, "genDimuon_genPtDR", dimuonPtBins, " #mu#mu pT (GeV)", dimuonDRBins, "#mu#mu #Delta R at IP");
245  myBook2D(ibooker, "globDimuon_genPtDR", dimuonPtBins, " #mu#mu pT (GeV)", dimuonDRBins, "#mu#mu #Delta R at IP");
246  myBook2D(ibooker, "globDimuon_recoPtDR", dimuonPtBins, " #mu#mu pT (GeV)", dimuonDRBins, "#mu#mu #Delta R at IP");
247  for(size_t i=0; i<filterNamesLevels.size(); i++){
248  myBook2D(ibooker, TString::Format("filt%dDimuon_recoPtDR",int(i+1)), dimuonPtBins, " #mu#mu pT (GeV)", dimuonDRBins, "#mu#mu #Delta R at IP", filterNamesLevels[i].first);
249  }
250  myBook2D(ibooker, "pathDimuon_recoPtDR", dimuonPtBins, " #mu#mu pT (GeV)", dimuonDRBins, "#mu#mu #Delta R at IP", triggerPathName);
251  for(size_t i=0; i<filterNamesLevels.size(); i++){
252  myBook2D(ibooker, TString::Format("diFilt%dDimuon_recoPtDR",int(i+1)), dimuonPtBins, " #mu#mu pT (GeV)", dimuonDRBins, "#mu#mu #Delta R at IP", filterNamesLevels[i].first);
253  }
254  myBook2D(ibooker, "diPathDimuon_recoPtDR", dimuonPtBins, " #mu#mu pT (GeV)", dimuonDRBins, "#mu#mu #Delta R at IP", triggerPathName);
255  myBook2D(ibooker, "resultDimuon_recoPtDR", dimuonPtBins, " #mu#mu pT (GeV)", dimuonDRBins, "#mu#mu #Delta R at IP");
256  // Pt DRpos Double
257  myBook2D(ibooker, "globDimuon_recoPtDRpos", dimuonPtBins, " #mu#mu pT (GeV)", dimuonDRBins, "#mu#mu #Delta R in MS");
258  for(size_t i=0; i<filterNamesLevels.size(); i++){
259  myBook2D(ibooker, TString::Format("filt%dDimuon_recoPtDRpos",int(i+1)), dimuonPtBins, " #mu#mu pT (GeV)", dimuonDRBins, "#mu#mu #Delta R in MS", filterNamesLevels[i].first);
260  }
261  myBook2D(ibooker, "pathDimuon_recoPtDRpos", dimuonPtBins, " #mu#mu pT (GeV)", dimuonDRBins, "#mu#mu #Delta R in MS", triggerPathName);
262  for(size_t i=0; i<filterNamesLevels.size(); i++){
263  myBook2D(ibooker, TString::Format("diFilt%dDimuon_recoPtDRpos",int(i+1)), dimuonPtBins, " #mu#mu pT (GeV)", dimuonDRBins, "#mu#mu #Delta R in MS", filterNamesLevels[i].first);
264  }
265  myBook2D(ibooker, "diPathDimuon_recoPtDRpos", dimuonPtBins, " #mu#mu pT (GeV)", dimuonDRBins, "#mu#mu #Delta R in MS", triggerPathName);
266  myBook2D(ibooker, "resultDimuon_recoPtDRpos", dimuonPtBins, " #mu#mu pT (GeV)", dimuonDRBins, "#mu#mu #Delta R in MS");
267 
268  // Matching
269  myBook2D(ibooker, "globGen_deltaEtaDeltaPhi", deltaEtaBins, "#Delta eta", deltaPhiBins, "#Delta phi");
270  for(size_t i=0; i<filterNamesLevels.size(); i++){
271  myBook2D(ibooker, TString::Format("filt%dGlob_deltaEtaDeltaPhi",int(i+1)), deltaEtaBins, "#Delta eta", deltaPhiBins, "#Delta phi", filterNamesLevels[i].first);
272  }
273  myBook2D(ibooker, "pathGlob_deltaEtaDeltaPhi", deltaEtaBins, "#Delta eta", deltaPhiBins, "#Delta phi", triggerPathName);
274  // Size of containers
275  vector<double> sizeBins; sizeBins.push_back(10); sizeBins.push_back(0); sizeBins.push_back(10);
276  myBook1D(ibooker, "genMuon_size", sizeBins, "container size" );
277  myBook1D(ibooker, "globMuon_size", sizeBins, "container size" );
278  for(size_t i=0; i<filterNamesLevels.size(); i++){
279  myBook1D(ibooker, TString::Format("filt%dMuon_size",int(i+1)), sizeBins, "container size", filterNamesLevels[i].first);
280  }
281  myBook1D(ibooker, "pathMuon_size", sizeBins, "container size", triggerPathName );
282 }
283 
285  if(filterNamesLevels.size()==0){
286  return;
287  }
288 //access the containers and create LeafCandidate copies
289  vector<LeafCandidate> genMuons;
291  iEvent.getByToken(genParticlesToken, genParticles);
292  if(genParticles.isValid()){
293  for(GenParticleCollection::const_iterator p=genParticles->begin(); p!= genParticles->end(); ++p){
294  if( p->status() == 1 && std::abs(p->pdgId())==13 &&
295  ( find( motherIDs.begin(), motherIDs.end(), -1 )!=motherIDs.end() || find( motherIDs.begin(), motherIDs.end(), getMotherId(&(*p)) )!=motherIDs.end() ) ){
296  genMuons.push_back( *p );
297  }
298  }
299  }else{
300  LogDebug("HLTriggerOfflineHeavyFlavor")<<"Could not access GenParticleCollection"<<endl;
301  }
302  sort(genMuons.begin(), genMuons.end(), GreaterByPt<LeafCandidate>());
303  ME["genMuon_size"]->Fill(genMuons.size());
304  LogDebug("HLTriggerOfflineHeavyFlavor")<<"GenParticleCollection from "<<genParticlesTag<<" has size: "<<genMuons.size()<<endl;
305 
306  vector<LeafCandidate> globMuons;
307  vector<LeafCandidate> globMuons_position;
308  Handle<MuonCollection> recoMuonsHandle;
309  iEvent.getByToken(recoMuonsToken, recoMuonsHandle);
310  if(recoMuonsHandle.isValid()){
311  for(MuonCollection::const_iterator p=recoMuonsHandle->begin(); p!= recoMuonsHandle->end(); ++p){
312  if(p->isGlobalMuon()){
313  globMuons.push_back(*p);
314  globMuons_position.push_back( LeafCandidate( p->charge(), math::XYZTLorentzVector(p->outerTrack()->innerPosition().x(), p->outerTrack()->innerPosition().y(), p->outerTrack()->innerPosition().z(), 0.) ) );
315  }
316  }
317  }else{
318  LogDebug("HLTriggerOfflineHeavyFlavor")<<"Could not access reco Muons"<<endl;
319  }
320  ME["globMuon_size"]->Fill(globMuons.size());
321  LogDebug("HLTriggerOfflineHeavyFlavor")<<"Global Muons from "<<recoMuonsTag<<" has size: "<<globMuons.size()<<endl;
322 
323 // access RAW trigger event
324  vector<vector<LeafCandidate> > muonsAtFilter;
325  vector<vector<LeafCandidate> > muonPositionsAtFilter;
326  for(size_t i=0; i<filterNamesLevels.size(); i++){
327  muonsAtFilter.push_back(vector<LeafCandidate>());
328  muonPositionsAtFilter.push_back(vector<LeafCandidate>());
329  }
330  Handle<TriggerEventWithRefs> rawTriggerEvent;
331  iEvent.getByToken( triggerSummaryRAWTag, rawTriggerEvent );
332  if( rawTriggerEvent.isValid() ){
333  for(size_t i=0; i<filterNamesLevels.size(); i++){
334  size_t index = rawTriggerEvent->filterIndex(InputTag( filterNamesLevels[i].first, "", triggerProcessName ));
335  if ( index < rawTriggerEvent->size() ){
336  if( filterNamesLevels[i].second==1 ){
337  vector<L1MuonParticleRef> l1Cands;
338  rawTriggerEvent->getObjects( index, TriggerL1Mu, l1Cands );
339  for(size_t j=0; j<l1Cands.size(); j++){
340  muonsAtFilter[i].push_back(*l1Cands[j]);
341  }
342  }else{
343  vector<RecoChargedCandidateRef> hltCands;
344  rawTriggerEvent->getObjects( index, TriggerMuon, hltCands );
345  for(size_t j=0; j<hltCands.size(); j++){
346  muonsAtFilter[i].push_back(*hltCands[j]);
347  if( filterNamesLevels[i].second==2 ){
348  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.) ) );
349  }
350  }
351  }
352  }
353  ME[TString::Format("filt%dMuon_size",int(i+1))]->Fill(muonsAtFilter[i].size());
354  LogDebug("HLTriggerOfflineHeavyFlavor")<<"Filter \""<<filterNamesLevels[i].first<<"\" has "<<muonsAtFilter[i].size()<<" muons"<<endl;
355  }
356  }else{
357  LogDebug("HLTriggerOfflineHeavyFlavor")<<"Could not access RAWTriggerEvent"<<endl;
358  }
359 
360 // access AOD trigger event
361  vector<LeafCandidate> pathMuons;
362  Handle<TriggerEvent> aodTriggerEvent;
363  iEvent.getByToken(triggerSummaryAODTag,aodTriggerEvent);
364  if(aodTriggerEvent.isValid()){
365  TriggerObjectCollection allObjects = aodTriggerEvent->getObjects();
366  for(int i=0; i<aodTriggerEvent->sizeFilters(); i++){
367  if(aodTriggerEvent->filterTag(i)==InputTag((filterNamesLevels.end()-1)->first,"",triggerProcessName)){
368  Keys keys = aodTriggerEvent->filterKeys(i);
369  for(size_t j=0; j<keys.size(); j++){
370  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 ) ) );
371  }
372  }
373  }
374  ME["pathMuon_size"]->Fill(pathMuons.size());
375  LogDebug("HLTriggerOfflineHeavyFlavor")<<"Path \""<<triggerPathName<<"\" has "<<pathMuons.size()<<" muons at last filter \""<<(filterNamesLevels.end()-1)->first<<"\""<<endl;
376  }else{
377  LogDebug("HLTriggerOfflineHeavyFlavor")<<"Could not access AODTriggerEvent"<<endl;
378  }
379 
380 // access Trigger Results
381  bool triggerFired = false;
383  iEvent.getByToken(triggerResultsToken,triggerResults);
384  if(triggerResults.isValid()){
385  LogDebug("HLTriggerOfflineHeavyFlavor")<<"Successfully initialized "<<triggerResultsTag<<endl;
386  const edm::TriggerNames & triggerNames = iEvent.triggerNames(*triggerResults);
387  bool hlt_exists = false;
388  for ( unsigned int i=0; i!=triggerNames.size(); i++) {
389  TString hlt_name = triggerNames.triggerName(i);
390  if (hlt_name.Contains(triggerPathName)) {
391  triggerFired = triggerResults->accept( i );
392  hlt_exists = true;
393  break;
394  }
395  }
396  if (!hlt_exists) {
397  LogDebug("HLTriggerOfflineHeavyFlavor")<<triggerResultsTag<<" has no trigger: "<<triggerPathName<<endl;
398  }
399  }else{
400  LogDebug("HLTriggerOfflineHeavyFlavor")<<"Could not initialize "<<triggerResultsTag<<endl;
401  }
402 
403 //create matching maps
404  vector<int> glob_gen(genMuons.size(),-1);
405  match( ME["globGen_deltaEtaDeltaPhi"], genMuons, globMuons, genGlobDeltaRMatchingCut, glob_gen );
406  vector<vector<int> > filt_glob;
407  for(size_t i=0; i<filterNamesLevels.size(); i++){
408  filt_glob.push_back( vector<int>(globMuons.size(),-1) );
409  if( filterNamesLevels[i].second == 1 ){
410  match( ME[TString::Format("filt%dGlob_deltaEtaDeltaPhi",int(i+1))], globMuons_position, muonsAtFilter[i] ,globL1DeltaRMatchingCut, filt_glob[i] );
411  }else if( filterNamesLevels[i].second == 2 ){
412  match( ME[TString::Format("filt%dGlob_deltaEtaDeltaPhi",int(i+1))], globMuons_position, muonPositionsAtFilter[i] ,globL2DeltaRMatchingCut, filt_glob[i] );
413  }else if( filterNamesLevels[i].second > 2){
414  match( ME[TString::Format("filt%dGlob_deltaEtaDeltaPhi",int(i+1))], globMuons, muonsAtFilter[i] ,globL3DeltaRMatchingCut, filt_glob[i] );
415  }
416  }
417  vector<int> path_glob(globMuons.size(),-1);
418  if( (filterNamesLevels.end()-1)->second == 1 ){
419  match( ME["pathGlob_deltaEtaDeltaPhi"], globMuons_position, pathMuons ,globL1DeltaRMatchingCut, path_glob );
420  }else if( (filterNamesLevels.end()-1)->second == 2 ){
421  match( ME["pathGlob_deltaEtaDeltaPhi"], globMuons, pathMuons ,globL2DeltaRMatchingCut, path_glob );
422  }else if( (filterNamesLevels.end()-1)->second > 2){
423  match( ME["pathGlob_deltaEtaDeltaPhi"], globMuons, pathMuons ,globL3DeltaRMatchingCut, path_glob );
424  }
425 
426 //fill histos
427  bool first = true;
428  for(size_t i=0; i<genMuons.size(); i++){
429  ME["genMuon_genEtaPt"]->Fill(genMuons[i].eta(), genMuons[i].pt());
430  ME["genMuon_genEtaPhi"]->Fill(genMuons[i].eta(), genMuons[i].phi());
431  if(glob_gen[i] != -1) {
432  ME["resGlobGen_genEtaPt"]->Fill(genMuons[i].eta(), genMuons[i].pt(), (globMuons[glob_gen[i]].pt()-genMuons[i].pt())/genMuons[i].pt() );
433  ME["globMuon_genEtaPt"]->Fill(genMuons[i].eta(), genMuons[i].pt());
434  ME["globMuon_genEtaPhi"]->Fill(genMuons[i].eta(), genMuons[i].phi());
435  ME["globMuon_recoEtaPt"]->Fill(globMuons[glob_gen[i]].eta(), globMuons[glob_gen[i]].pt());
436  ME["globMuon_recoEtaPhi"]->Fill(globMuons[glob_gen[i]].eta(), globMuons[glob_gen[i]].phi());
437  for(size_t f=0; f<filterNamesLevels.size(); f++) {
438  if(filt_glob[f][glob_gen[i]] != -1) {
439  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() );
440  ME[TString::Format("filt%dMuon_recoEtaPt",int(f+1))]->Fill(globMuons[glob_gen[i]].eta(), globMuons[glob_gen[i]].pt());
441  ME[TString::Format("filt%dMuon_recoEtaPhi",int(f+1))]->Fill(globMuons[glob_gen[i]].eta(), globMuons[glob_gen[i]].phi());
442  }else{
443  break;
444  }
445  }
446  if(path_glob[glob_gen[i]] != -1){
447  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() );
448  ME["pathMuon_recoEtaPt"]->Fill(globMuons[glob_gen[i]].eta(), globMuons[glob_gen[i]].pt());
449  ME["pathMuon_recoEtaPhi"]->Fill(globMuons[glob_gen[i]].eta(), globMuons[glob_gen[i]].phi());
450  }
451 //highest pt muon
452  if( first ){
453  first = false;
454  if( triggerFired ){
455  ME["resultMuon_recoEtaPt"]->Fill(globMuons[glob_gen[i]].eta(), globMuons[glob_gen[i]].pt());
456  ME["resultMuon_recoEtaPhi"]->Fill(globMuons[glob_gen[i]].eta(), globMuons[glob_gen[i]].phi());
457  }
458  }
459  }
460  }
461 
462 //fill dimuon histograms (highest pT, opposite charge)
463  int secondMuon = 0;
464  for(size_t j=1; j<genMuons.size(); j++){
465  if(genMuons[0].charge()*genMuons[j].charge()==-1){
466  secondMuon = j;
467  break;
468  }
469  }
470  if(secondMuon > 0){
471 //two generated
472  double genDimuonPt = (genMuons[0].p4()+genMuons[secondMuon].p4()).pt();
473  double genDimuonEta = (genMuons[0].p4()+genMuons[secondMuon].p4()).eta();
474  double genDimuonRap = (genMuons[0].p4()+genMuons[secondMuon].p4()).Rapidity();
475  double genDimuonDR = deltaR<LeafCandidate,LeafCandidate>( genMuons[0], genMuons[secondMuon] );
476  bool highPt = genMuons[0].pt()>7. && genMuons[secondMuon].pt()>7;
477  ME["genDimuon_genEtaPt"]->Fill( genDimuonEta, genDimuonPt );
478  ME["genDimuon_genRapPt"]->Fill( genDimuonRap, genDimuonPt );
479  if(highPt) ME["genDimuon_genPtDR"]->Fill( genDimuonPt, genDimuonDR );
480 //two global
481  if(glob_gen[0]!=-1 && glob_gen[secondMuon]!=-1){
482  ME["globDimuon_genEtaPt"]->Fill( genDimuonEta, genDimuonPt );
483  ME["globDimuon_genRapPt"]->Fill( genDimuonRap, genDimuonPt );
484  if(highPt) ME["globDimuon_genPtDR"]->Fill( genDimuonPt, genDimuonDR );
485  double globDimuonPt = (globMuons[glob_gen[0]].p4()+globMuons[glob_gen[secondMuon]].p4()).pt();
486  double globDimuonEta = (globMuons[glob_gen[0]].p4()+globMuons[glob_gen[secondMuon]].p4()).eta();
487  double globDimuonRap = (globMuons[glob_gen[0]].p4()+globMuons[glob_gen[secondMuon]].p4()).Rapidity();
488  double globDimuonDR = deltaR<LeafCandidate,LeafCandidate>( globMuons[glob_gen[0]], globMuons[glob_gen[secondMuon]] );
489  double globDimuonDRpos = deltaR<LeafCandidate,LeafCandidate>( globMuons_position[glob_gen[0]], globMuons_position[glob_gen[secondMuon]] );
490  ME["globDimuon_recoEtaPt"]->Fill( globDimuonEta, globDimuonPt );
491  ME["globDimuon_recoRapPt"]->Fill( globDimuonRap, globDimuonPt );
492  if(highPt) ME["globDimuon_recoPtDR"]->Fill( globDimuonPt, globDimuonDR );
493  if(highPt) ME["globDimuon_recoPtDRpos"]->Fill( globDimuonPt, globDimuonDRpos );
494 //two filter objects
495  for(size_t f=0; f<filterNamesLevels.size(); f++){
496  if(filt_glob[f][glob_gen[0]] != -1 && filt_glob[f][glob_gen[secondMuon]] != -1){
497  ME[TString::Format("diFilt%dDimuon_recoEtaPt",int(f+1))]->Fill( globDimuonEta, globDimuonPt );
498  ME[TString::Format("diFilt%dDimuon_recoRapPt",int(f+1))]->Fill( globDimuonRap, globDimuonPt );
499  if(highPt) ME[TString::Format("diFilt%dDimuon_recoPtDR",int(f+1))]->Fill( globDimuonPt, globDimuonDR );
500  if(highPt) ME[TString::Format("diFilt%dDimuon_recoPtDRpos",int(f+1))]->Fill( globDimuonPt, globDimuonDRpos );
501  }else{
502  break;
503  }
504  }
505 //one filter object
506  for(size_t f=0; f<filterNamesLevels.size(); f++){
507  if(filt_glob[f][glob_gen[0]] != -1 || filt_glob[f][glob_gen[secondMuon]] != -1){
508  ME[TString::Format("filt%dDimuon_recoEtaPt",int(f+1))]->Fill( globDimuonEta, globDimuonPt );
509  ME[TString::Format("filt%dDimuon_recoRapPt",int(f+1))]->Fill( globDimuonRap, globDimuonPt );
510  if(highPt) ME[TString::Format("filt%dDimuon_recoPtDR",int(f+1))]->Fill( globDimuonPt, globDimuonDR );
511  if(highPt) ME[TString::Format("filt%dDimuon_recoPtDRpos",int(f+1))]->Fill( globDimuonPt, globDimuonDRpos );
512  }else{
513  break;
514  }
515  }
516 //two path objects
517  if(path_glob[glob_gen[0]] != -1 && path_glob[glob_gen[secondMuon]] != -1){
518  ME["diPathDimuon_recoEtaPt"]->Fill( globDimuonEta, globDimuonPt );
519  ME["diPathDimuon_recoRapPt"]->Fill( globDimuonRap, globDimuonPt );
520  if(highPt) ME["diPathDimuon_recoPtDR"]->Fill( globDimuonPt, globDimuonDR );
521  if(highPt) ME["diPathDimuon_recoPtDRpos"]->Fill( globDimuonPt, globDimuonDRpos );
522  }
523 //one path object
524  if(path_glob[glob_gen[0]] != -1 || path_glob[glob_gen[secondMuon]] != -1){
525  ME["pathDimuon_recoEtaPt"]->Fill( globDimuonEta, globDimuonPt );
526  ME["pathDimuon_recoRapPt"]->Fill( globDimuonRap, globDimuonPt );
527  if(highPt) ME["pathDimuon_recoPtDR"]->Fill( globDimuonPt, globDimuonDR );
528  if(highPt) ME["pathDimuon_recoPtDRpos"]->Fill( globDimuonPt, globDimuonDRpos );
529  }
530 //trigger result
531  if( triggerFired ){
532  ME["resultDimuon_recoEtaPt"]->Fill( globDimuonEta, globDimuonPt );
533  ME["resultDimuon_recoRapPt"]->Fill( globDimuonRap, globDimuonPt );
534  if(highPt) ME["resultDimuon_recoPtDR"]->Fill( globDimuonPt, globDimuonDR );
535  if(highPt) ME["resultDimuon_recoPtDRpos"]->Fill( globDimuonPt, globDimuonDRpos );
536  }
537  }
538  }
539 }
540 
542  const Candidate* mother = p->mother();
543  if( mother ){
544  if( mother->pdgId() == p->pdgId() ){
545  return getMotherId(mother);
546  }else{
547  return mother->pdgId();
548  }
549  }else{
550  return 0;
551  }
552 }
553 
554 void HeavyFlavorValidation::match( MonitorElement * me, vector<LeafCandidate> & from, vector<LeafCandidate> & to, double dRMatchingCut, vector<int> & map ){
555  vector<double> dR(from.size());
556  for(size_t i=0; i<from.size(); i++){
557  map[i] = -1;
558  dR[i] = 10.;
559  //find closest
560  for(size_t j=0; j<to.size(); j++){
561  double dRtmp = deltaR<double>(from[i].eta(), from[i].phi(), to[j].eta(), to[j].phi());
562  if( dRtmp < dR[i] ){
563  dR[i] = dRtmp;
564  map[i] = j;
565  }
566  }
567  //fill matching histo
568  if( map[i] != -1 ){
569  me->Fill( to[map[i]].eta()-from[i].eta(), deltaPhi<double>(to[map[i]].phi(), from[i].phi()) );
570  }
571  //apply matching cut
572  if( dR[i] > dRMatchingCut ){
573  map[i] = -1;
574  }
575  //remove duplication
576  if( map[i] != -1 ){
577  for(size_t k=0; k<i; k++){
578  if( map[k] != -1 && map[i] == map[k] ){
579  if( dR[i] < dR[k] ){
580  map[k] = -1;
581  }else{
582  map[i] = -1;
583  }
584  break;
585  }
586  }
587  }
588  }
589 }
590 
591 void HeavyFlavorValidation::myBook2D(DQMStore::IBooker & ibooker, TString name, vector<double> &ptBins, TString ptLabel, vector<double> &etaBins, TString etaLabel, TString title )
592 {
593 // dqmStore->setCurrentFolder(dqmFolder+"/"+folder);
594  int ptN = ptBins.size()==3 ? (int)ptBins[0]+1 : ptBins.size();
595  Double_t *pt = new Double_t[ ptN ];
596  for(int i=0; i<ptN; i++){
597  pt[i] = ptBins.size()==3 ? ptBins[1] + i*(ptBins[2]-ptBins[1])/ptBins[0] : ptBins[i] ;
598  }
599  int etaN = etaBins.size()==3 ? (int)etaBins[0]+1 : etaBins.size();
600  Double_t *eta = new Double_t[ etaN ];
601  for(int i=0; i<etaN; i++){
602  eta[i] = etaBins.size()==3 ? etaBins[1] + i*(etaBins[2]-etaBins[1])/etaBins[0] : etaBins[i] ;
603  }
604  TH2F *h = new TH2F( name, name, ptN-1, pt, etaN-1, eta );
605  h->SetXTitle(ptLabel);
606  h->SetYTitle(etaLabel);
607  h->SetTitle(title);
608  ME[name] = ibooker.book2D( name.Data(), h );
609  delete h;
610 }
611 
612 void HeavyFlavorValidation::myBookProfile2D(DQMStore::IBooker & ibooker, TString name, vector<double> &ptBins, TString ptLabel, vector<double> &etaBins, TString etaLabel, TString title )
613 {
614 // dqmStore->setCurrentFolder(dqmFolder+"/"+folder);
615  int ptN = ptBins.size()==3 ? (int)ptBins[0]+1 : ptBins.size();
616  Double_t *pt = new Double_t[ ptN ];
617  for(int i=0; i<ptN; i++){
618  pt[i] = ptBins.size()==3 ? ptBins[1] + i*(ptBins[2]-ptBins[1])/ptBins[0] : ptBins[i] ;
619  }
620  int etaN = etaBins.size()==3 ? (int)etaBins[0]+1 : etaBins.size();
621  Double_t *eta = new Double_t[ etaN ];
622  for(int i=0; i<etaN; i++){
623  eta[i] = etaBins.size()==3 ? etaBins[1] + i*(etaBins[2]-etaBins[1])/etaBins[0] : etaBins[i] ;
624  }
625  TProfile2D *h = new TProfile2D( name, name, ptN-1, pt, etaN-1, eta );
626  h->SetXTitle(ptLabel);
627  h->SetYTitle(etaLabel);
628  h->SetTitle(title);
629  ME[name] = ibooker.bookProfile2D( name.Data(), h );
630  delete h;
631 }
632 
633 void HeavyFlavorValidation::myBook1D(DQMStore::IBooker & ibooker, TString name, vector<double> &bins, TString label, TString title )
634 {
635 // dqmStore->setCurrentFolder(dqmFolder+"/"+folder);
636  int binsN = bins.size()==3 ? (int)bins[0]+1 : bins.size();
637  Double_t *myBins = new Double_t[ binsN ];
638  for(int i=0; i<binsN; i++){
639  myBins[i] = bins.size()==3 ? bins[1] + i*(bins[2]-bins[1])/bins[0] : bins[i] ;
640  }
641  TH1F *h = new TH1F( name, name, binsN-1, myBins );
642  h->SetXTitle(label);
643  h->SetTitle(title);
644  ME[name] = ibooker.book1D( name.Data(), h );
645  delete h;
646 }
647 
649 }
650 
651 //define this as a plug-in
#define LogDebug(id)
T getParameter(std::string const &) const
T getUntrackedParameter(std::string const &, T const &) const
int i
Definition: DBlmapReader.cc:9
virtual edm::TriggerNames const & triggerNames(edm::TriggerResults const &triggerResults) const
Definition: Event.cc:220
void myBook2D(DQMStore::IBooker &ibooker, TString name, vector< double > &xBins, TString xLabel, vector< double > &yBins, TString yLabel, TString title)
virtual const Candidate * mother(size_type i=0) const =0
return pointer to mother
virtual void analyze(const edm::Event &, const edm::EventSetup &) override
bool getByToken(EDGetToken token, Handle< PROD > &result) const
Definition: Event.h:464
enum start value shifted to 81 so as to avoid clashes with PDG codes
void bookHistograms(DQMStore::IBooker &, edm::Run const &, edm::EventSetup const &) override
void cd(void)
Definition: DQMStore.cc:266
#define DEFINE_FWK_MODULE(type)
Definition: MakerMacros.h:17
const std::vector< std::string > & triggerNames() const
names of trigger paths
vector< pair< string, int > > filterNamesLevels
const double muonMass
Strings::size_type size() const
Definition: TriggerNames.cc:39
void find(edm::Handle< EcalRecHitCollection > &hits, DetId thisDet, std::vector< EcalRecHitCollection::const_iterator > &hit, bool debug=false)
Definition: FindCaloHit.cc:7
void bookHistograms(fwlite::EventContainer &eventCont)
EDGetTokenT< TriggerEvent > triggerSummaryAODTag
void Fill(long long x)
U second(std::pair< T, U > const &p)
PtEtaPhiMLorentzVectorD PtEtaPhiMLorentzVector
Lorentz vector with cartesian internal representation.
Definition: LorentzVector.h:25
XYZTLorentzVectorD XYZTLorentzVector
Lorentz vector with cylindrical internal representation using pseudorapidity.
Definition: LorentzVector.h:29
Definition: ME.h:11
void myBookProfile2D(DQMStore::IBooker &ibooker, TString name, vector< double > &xBins, TString xLabel, vector< double > &yBins, TString yLabel, TString title)
int getMotherId(const Candidate *p)
int iEvent
Definition: GenABIO.cc:230
void match(MonitorElement *me, vector< LeafCandidate > &from, vector< LeafCandidate > &to, double deltaRMatchingCut, vector< int > &map)
MonitorElement * bookProfile2D(Args &&...args)
Definition: DQMStore.h:163
EDGetTokenT< TriggerEventWithRefs > triggerSummaryRAWTag
MonitorElement * book1D(Args &&...args)
Definition: DQMStore.h:115
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
int j
Definition: DBlmapReader.cc:9
void dqmBeginRun(const edm::Run &, const edm::EventSetup &)
double f[11][100]
void myBook1D(DQMStore::IBooker &ibooker, TString name, vector< double > &xBins, TString label, TString title)
The Signals That Services Can Subscribe To This is based on ActivityRegistry h
Helper function to determine trigger accepts.
Definition: Activities.doc:4
static std::string const triggerResults
Definition: EdmProvDump.cc:40
bool isValid() const
Definition: HandleBase.h:75
HeavyFlavorValidation(const edm::ParameterSet &)
const std::vector< std::string > & moduleLabels(unsigned int trigger) const
label(s) of module(s) on a trigger path
virtual int pdgId() const =0
PDG identifier.
std::vector< TriggerObject > TriggerObjectCollection
collection of trigger physics objects (e.g., all isolated muons)
Definition: TriggerObject.h:81
void setCurrentFolder(const std::string &fullpath)
Definition: DQMStore.cc:274
MonitorElement * book2D(Args &&...args)
Definition: DQMStore.h:133
void myBook2D(DQMStore::IBooker &ibooker, TString name, vector< double > &xBins, TString xLabel, vector< double > &yBins, TString yLabel)
EDGetTokenT< MuonCollection > recoMuonsToken
std::string const & triggerName(unsigned int index) const
Definition: TriggerNames.cc:27
std::vector< size_type > Keys
bool init(const edm::Run &iRun, const edm::EventSetup &iSetup, const std::string &processName, bool &changed)
d&#39;tor
EDGetTokenT< GenParticleCollection > genParticlesToken
void myBookProfile2D(DQMStore::IBooker &ibooker, TString name, vector< double > &xBins, TString xLabel, vector< double > &yBins, TString yLabel)
std::pair< typename Association::data_type::first_type, double > match(Reference key, Association association, bool bestMatchByMaxValue)
Generic matching function.
Definition: Utils.h:10
tuple level
Definition: testEve_cfg.py:34
void myBook1D(DQMStore::IBooker &ibooker, TString name, vector< double > &xBins, TString label)
EDGetTokenT< TriggerResults > triggerResultsToken
map< TString, MonitorElement * > ME
tuple size
Write out results.
Definition: Run.h:43